Flexible vocal timing adaptation to predictable and unpredictable playback noise

About

In this experiment, we played amplitude modulated white noise to a group of bats in order to investigate their ability to temporally modulate their calling behavior at the group level.

We hypothesized that bats would adapt the timing of their calls to avoid temporal overlap with the amplitude modulated noise.

After extracting call timins from acoustic recordings (.wav files) and identifying the instantaneous period (or phase) in the real or simulated (as in the silent condition) amplitude modulation cycle, we analyzed the call onset data, primarily with circular statistics.

Setup

packages <- c("circular", "CircStats", "CircSpaceTime", "sjPlot", "yarrr", "janitor", "broom", "knitr",  "kableExtra", "car", "MASS", "emmeans", "rcompanion", "pscl", "randomForest", "chisq.posthoc.test", "viridis", "scales","ggh4x","ggside", "ggdist", "ggridges","gganimate","plotly","htmlwidgets","tidyverse")

# Install packages not yet installed
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
  install.packages(packages[!installed_packages])
}

# Packages loading
invisible(lapply(packages, library, character.only = TRUE))

Parameters

str(params)

## List of 7
##  $ subsample   : logi FALSE
##  $ data_path   : chr "./../1-data"
##  $ results_path: chr "./../3-results"
##  $ exp         : chr "exp_2"
##  $ plot_sm     : logi TRUE
##  $ use_all     : logi TRUE
##  $ seed        : num 42

# local paths
data_path <- params$data_path
results_path <- params$results_path
tables_path <- file.path(results_path,'tables')
  if (!dir.exists(tables_path)) dir.create(tables_path)
figures_path <- file.path(results_path,'figures')
  if (!dir.exists(figures_path)) dir.create(figures_path)
models_path <- file.path(results_path,'models')
  if (!dir.exists(models_path)) dir.create(models_path)
sm_path <- file.path(results_path,'supplementary_figures')
  if (!dir.exists(sm_path)) dir.create(sm_path)

# experiment number
exp <- params$exp

# seed for subsampling
seed <- params$seed

# whether to use all data for circular stats
use_all <- params$use_all
subsample <- params$subsample

# whether to plot supplementary plots/tables
supplementary <- params$plot_sm

Load data

# load data
data <- readr::read_rds(file.path(data_path,'input',paste0(exp, "_data.RDS"))) %>%
  mutate(start_phase_circ = circular(start_phase, units="radians", zero = 0, modulo = "2pi", rotation = "counter"),.after ="start_phase") 

# data for averaging histograms (modulation variable is numeric-ish)
data_2proc <- data %>%
  mutate(modulation = factor(modulation, 
                             levels = sort(unique(data$modulation))))

# data for remaining analyses (modulation variable is character-ish, better for plotting)
data <- data %>%
  mutate(modulation = factor(modulation, 
                             levels = sort(unique(data$modulation)), 
                             labels = paste0(as.character(sort(unique(data$modulation))),"Hz")))
str(data)

## 'data.frame':    1088994 obs. of  23 variables:
##  $ group           : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
##  $ session         : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ part            : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ condition       : Factor w/ 3 levels "silence","steady-state masker",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ modulation      : Factor w/ 8 levels "4Hz","8Hz","16Hz",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ minute          : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ start_seconds   : num  0.0151 0.0777 0.1069 0.1599 0.1902 ...
##  $ stop_seconds    : num  0.0169 0.0794 0.1086 0.1627 0.1918 ...
##  $ start_sp        : num  3776 19437 26717 39986 47554 ...
##  $ end_sp          : num  4223 19848 27141 40671 47960 ...
##  $ start_cycle     : int  1 2 2 3 4 4 5 6 7 9 ...
##  $ end_cycle       : int  1 2 2 3 4 4 5 6 7 9 ...
##  $ start_phase     : num  1.518 1.533 4.46 3.513 0.273 ...
##  $ start_phase_circ: 'circular' num  1.518 1.533 4.46 3.513 0.273 ...
##   ..- attr(*, "circularp")=List of 6
##   .. ..$ type    : chr "angles"
##   .. ..$ units   : chr "radians"
##   .. ..$ template: chr "none"
##   .. ..$ modulo  : chr "2pi"
##   .. ..$ zero    : num 0
##   .. ..$ rotation: chr "counter"
##  $ end_phase       : num  1.698 1.698 4.631 3.788 0.436 ...
##  $ start_per       : num  0.0151 0.01524 0.04437 0.03494 0.00271 ...
##  $ end_per         : num  0.01689 0.01689 0.04606 0.03768 0.00434 ...
##  $ duration        : num  0.00179 0.00164 0.0017 0.00274 0.00162 ...
##  $ onset_interval  : num  NA 0.0626 0.0291 0.0531 0.0303 ...
##  $ class           : chr  "call" "call" "call" "call" ...
##  $ file_no         : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ familiarization : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ order           : num  1 1 1 1 1 1 1 1 1 1 ...

# save interim data manips to 1-data/output
data_path <- file.path(data_path,'output')
  if (!dir.exists(data_path)) dir.create(data_path)

Summary Stats

A few rough summary statistics. Note that each row is an observation (call event), with parameters:

• group: experimental group (1:4) (Note that for experiment 1, there were 8 groups of bats tested in each modulation rate context, but each is labelled 1-4.)
• session: recording day (1:5)
• part: experimental block (1:3), silent condition is always part 1
• condition: playback condition
• modulation: amplitude modulation rate for masking conditions
• minute: minute since the start of the recording file in which the call occurred
• start_seconds: timestamps of call onsets in seconds since the start of the file
• stop_seconds: ibid. for call offsets
• start_sp: samples since start of file at call onset
• end_sp: ibid. for call offsets
• start_cycle: nth cycle of the corresponding modulation rate in the corresponding file at call onset
• end_cycle: ibid. at call onset
• start_phase: moment in the corresponding modulation cycle (in radians) at call onset
• start_phase_circ: ibid., but of class circular
• end_phase: ibid. at call offset
• start_per: timepoint in the corresponding modulation period (in seconds) at call onset
• end_per: ibid. at call offset
• duration: call offset - call onset (in seconds)
• onset_interval: call onset - previous call’s call onset (is NA if call is the first in the file)
• class: type of call (for now, all are labelled “call”)
• file_no: index of recording file, indicating the order of recorded files from the corresponding experimental block
• familiarization: number of days in which bats were housed together as a group prior to the first recording day
• order: randomization/counter-balancing order for the presentation of masking conditions

#summary(data)

Total number of calls:

nrow(data)

## [1] 1088994

#1,425,067 # exp 1
#1,088,994 # exp 2
# = 2,514,061 # total

Call durations:

# durations
# data1 <- data
# data2 <- data
# data <- rbind(data1,data2)
# nrow(data)

summary(data$duration)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000996 0.001716 0.003128 0.003881 0.004364 0.704744

# #    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
# #0.000996 0.001704 0.003372 0.004246 0.005036 0.704744 

IQR(data$duration)

## [1] 0.002648

# # 0.003332

Distribution of call onsets in the modulation cycle

Calculate histogram averages (essentially, a PSTH)

call_summary1 <- get_avg_hist(data_2proc, NULL, list("group","session"), "period", 
                              by_group = TRUE, by_session = TRUE)
call_summary60 <- get_avg_hist(data_2proc, 60, list("group","session"), "period")

# optionally: check that the two have the correct # of bins
# call_summary1$generic %>% group_by(modulation) %>% summarise(sort(unique(mids),decreasing = T)[2])
# call_summary60$generic %>% group_by(modulation) %>% summarise(n=length(unique(mids)))

First, we asked: Did the presence of playback noise affect the distribution of emitted calls within the modulation cycle?

Here we plot an averaged histogram of detected call onsets (binned calls (1 ms bins) summed over all cycles, averaged over groups and recording days):

Points are mean # calls, shaded areas are +- sem…

# average call counts all at once, 1 ms bins
p <- p.hist_avg_dot_kwargs(dat = call_summary1$generic, facets = c("condition","modulation"))
 print(p)

 ggsave(file.path(figures_path,paste0(exp, '_fig1a.png') ), p, width = w, height = 5)

# average call counts with fixed y-axis, to better illustrate the relative calling rate between conditions
p2 <- p.hist_avg_dot_kwargs_fixedy(dat = call_summary1$generic, facets = c("","modulation")) 
  print(p2)

 ggsave(file.path(figures_path,paste0(exp, '_fig1a-alt1.png') ), p2, width = w, height = 5)

# average call counts with fixed y-axis, with 60 ms bins for each modulation rate, to show trend on the same x-scale across conditions  
p3 <- p.hist_avg_dot_kwargs_fixedy(dat = call_summary60$generic, facets = c("","modulation")) 
  print(p3)

  ggsave(file.path(figures_path,paste0(exp, '_fig1a-alt2.png') ), p3, width = w, height = 5) 

# show average call counts (with and without fixed y axis) for each group, and for each recording day ("session")
for (m in sort(unique(call_summary1$group$modulation))) {
 p1 <- p.hist_avg_dot_kwargs(dat = call_summary1$group %>% filter(modulation == m), 
                             facets = c("condition","group")) +
                             ggtitle(paste0(m, ": by group"))
 p2 <- p.hist_avg_dot_kwargs_fixedy(dat = call_summary1$group %>% filter(modulation == m), 
                                    facets = c("","group")) +
                             ggtitle(paste0(m, ": by group")) 
 p3 <- p.hist_avg_dot_kwargs(dat = call_summary1$session %>% filter(modulation == m),
                             facets = c("condition","session")) +
                             ggtitle(paste0(m, ": by day"))
 
 print(p1); print(p2); print(p3)
 ggsave(file.path(sm_path,paste0(exp, '_figs3a_',m,'.png') ), 
        p1, width = 7, height = 5)
 ggsave(file.path(sm_path,paste0(exp, '_figs3a_',m,'-alt1.png') ), 
        p2, width = 7, height = 5)
 ggsave(file.path(sm_path,paste0(exp, '_figs3b_',m,'.png') ), 
        p3, width = 7, height = 5)
}

# show average hist of binned calls for the first 7.5-8 mins of continuous AM playback for each mod.
call_summary1_init <- get_avg_hist(data_init, NULL, list("group"), "period", by_group = TRUE)
#str(data_init)
 
firsts <- p.hist_avg_col(dat = call_summary1_init$group %>% 
                 dplyr::filter(condition=="full-band masker" | condition=="steady-state masker"), 
               facets = c("modulation","group"))

print(firsts)

ggsave(file.path(sm_path,paste0(exp, '_figs4.png') ), firsts, width = 7, height = h2)

Circular Statistics

To analyze the distribution of call onsets within the modulation cycle, we must use circular statistics (with functions from the package circular), since the nature of our data is cyclical.

First we will compute the following summary statistics:

• Mean Resultant Direction (\(\bar{\theta}\))
• Mean Resultant Length (\(\bar{R}\))
• Circular Variance (\(Vm\)): \(1-R\)
• Circular Standard Deviation (\(v\)): \(\sqrt{1-2*log(\bar{R})}\)
  
• MLE Mu (\(\mu\)): mean parameter estimated from the maximum likelihood von Mises distribution (circular::mle.vonmises())
• MLE Kappa (\(\kappa\)): concentration parameter estimated from the maximum likelihood von Mises distribution (circular::mle.vonmises())
• MLE CIs for \(\bar{\theta}\): confidence intervals for the mean parameters estimated from the maximum likelihood von Mises distribution (circular::mle.vonmises.bootstrap.ci())

Then, we’ll get bootsrapped mean and concentration parameters from the whole dataset.

---

We predicted that calls emitted in presence of modulated noise would tend to “cluster” unimodally in the hemisphere corresponding to the amplitude trough. Conversely, we predicted that the distribution of call onsets emitted in silence would show no clustering, and therefore resemble a uniform distribution.

This prediction entails two specific hypotheses:

1. Distribution of calls detected in silence will be equivalent to a uniform distribution, while that of calls detected in the presence of masking noise will deviate from the uniform distribution. –> Method: Rayleigh test of uniformity (circular::rayleigh.test(x, mu=circular(0)))

2. The distribution of calls in the playback conditions will each besignificantly different from the silence condition. –> Method: Mardia-Watson-Wheeler non-parametric test of homogeneous samples (circular::watson.wheeler.test(x,y)), where differences between samples can be in either the mean or the variance.

Calculate summary and distribution statistics

In which conditions were call onset distributions unimodally clustered?

# prepare data for analysis: split into lists, one for each modulation x condition case
# data is subsampled if requested in header

data_split <-  data_samp %>% group_by(modulation, condition) %>% group_split()

# label and convert to circular data
circ_data <- lapply(data_split, function(x) 
  list(modulation = unique(x$modulation),
       condition = unique(x$condition),
       phase = x$start_phase_circ)) 


if (file.exists(f_circsum)) { # as RDS...
  circ_summary <- readr::read_rds(f_circsum)
  
} else {
  # calculate values:
  circ_datsum <- lapply(circ_data, 
         function(x) tibble(modulation = x$modulation, 
                            condition = x$condition,
                            n = length(x$phase),
                            theta_bar = mean.circular(x$phase, na.rm = T), # angular mean
                            r_bar = rho.circular(x$phase, na.rm = T), # mean resultant length
                            vm = var.circular(x$phase), 
                            v = sd.circular(x$phase),  
                            # Test of uniformity:
                            rt = rayleigh.test(x$phase, mu=circular(0))$statistic,
                            p = rayleigh.test(x$phase, mu=circular(0))$p.value,
                            # MLE von Mises:
                            mle_mu = mle.vonmises(x$phase)$mu,
                            bs_mu_ci_l = mle.vonmises.bootstrap.ci(x$phase, alpha=.05)$mu.ci[1], 
                            bs_mu_ci_h = mle.vonmises.bootstrap.ci(x$phase, alpha=.05)$mu.ci[2],
                            mle_kappa = mle.vonmises(x$phase)$kappa) 
                        )  
  
  # wrangle & correct p values
  circ_summary <- tibble(do.call(rbind,circ_datsum))  %>%
                  mutate(p_adjust = p.adjust(p, method ="bonferroni"), .after = "p") %>%
                  mutate(across(where(is.numeric), round, 3))
  
  head(circ_summary)
  
  # save, if first time
  if (subsample) {
    write.csv(circ_summary, file = file.path(data_path,paste0(exp, "_circ_data_seed",seed,".csv")), row.names = F)
    readr::write_rds(circ_summary, file = file.path(data_path,paste0(exp, "_circ_data_seed",seed,".RDS")))
  } else {
    write.csv(circ_summary, file = file.path(data_path,paste0(exp, "_circ_data-alldata.csv")), row.names = F)
    readr::write_rds(circ_summary, file = file.path(data_path,paste0(exp, "_circ_data-alldata.RDS")))
  }

}

circ_summary

## # A tibble: 24 × 14
##    modulation condi…¹      n theta…² r_bar    vm     v    rt     p p_adj…³ mle_mu bs_mu…⁴ bs_mu…⁵ mle_k…⁶
##    <fct>      <fct>    <dbl> <circu> <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl> <circ> <circu> <circu>   <dbl>
##  1 4Hz        silence  59635 6.203   0.003 0.997  3.39 0.003 0.134       1 6.203  4.207   1.704     0.006
##  2 4Hz        steady…  28938 5.246   0.074 0.926  2.28 0.038 0           0 5.246  5.132   5.353     0.148
##  3 4Hz        random… 105990 5.245   0.054 0.946  2.42 0.027 0           0 5.245  5.165   5.326     0.108
##  4 8Hz        silence  62203 5.480   0.002 0.998  3.46 0.002 0.27        1 5.480  2.943   2.224     0.005
##  5 8Hz        steady…  65446 5.795   0.057 0.943  2.40 0.05  0           0 5.795  5.699   5.890     0.114
##  6 8Hz        random…  51939 5.459   0.061 0.939  2.37 0.041 0           0 5.459  5.366   5.567     0.122
##  7 16Hz       silence  65769 0.858   0.003 0.997  3.38 0.002 0.213       1 0.858  5.194   3.164     0.007
##  8 16Hz       steady…  60482 5.732   0.042 0.958  2.52 0.036 0           0 5.732  5.610   5.851     0.084
##  9 16Hz       random…  25543 5.376   0.046 0.954  2.49 0.028 0           0 5.376  5.178   5.571     0.091
## 10 25Hz       silence  67403 4.819   0.003 0.997  3.37 0     0.448       1 4.819  2.495   1.907     0.007
## # … with 14 more rows, and abbreviated variable names ¹​condition, ²​theta_bar, ³​p_adjust, ⁴​bs_mu_ci_l,
## #   ⁵​bs_mu_ci_h, ⁶​mle_kappa

Plot circular call density

Now, let’s plot our call count data (as a density rather than raw call count) as before, but on the polar plane and indicate the mean direction and resultant length of the resultant vector.

circ_density <- p.circ_density(dat = data_samp, circ_dat = circ_summary_plot, bins = 30) %>% recolor(which="f")

circ_density

if (subsample) {
  ggsave(file.path(figures_path,paste0(exp, "_fig1b_",seed,".png")), circ_density, width = w, height = 5)
} else {
  ggsave(file.path(figures_path,paste0(exp, '_fig1b.png')), circ_density, width = w, height = 5)
}

Calculate bootstrapped angular vectors

How consistent was the clustering, if present?

Plot bs parameters in the polar plane

circ_params <- p.circ_params(mu_bs_df, compress = T)  %>% recolor()


circ_params

ggsave(file.path(figures_path,paste0(exp, '_fig1c.png')), circ_params, width = 3, height = 6)

… and in the cartesian plane

mu_kappa <- p.mu_kappa(dat = mu_bs_df %>% 
                         group_by(modulation, condition) %>% 
                         slice_sample(n = 500)) %>% recolor()
mu_kappa

ggsave(file.path(figures_path,paste0(exp, '_fig1d-500.png') ), mu_kappa, width = 6, height = 6)

[Mardia]-Watson-Wheeler Homogeneity of Samples

We now compute frequentist non-parametric analyses for whether the distributions of calls observed in our different modulation and masking conditions are identical.

The null hypothesis is that the distributions are identical, either in their means or their variances:

Compare playback conditions within modulation rate contexts

Did distributions of call onsets differ between playback conditions at each modulation rate?

circ_data_df <- do.call(rbind, lapply(circ_data, function(x) as.data.frame(x)))

ww_tests <- lapply(sort(unique(circ_data_df$modulation)), function(x) 
  watson.wheeler.test(phase ~ condition,circ_data_df[circ_data_df$modulation==x,]) %>%
    broom::tidy() %>%
     mutate(across(where(is.numeric), round, 3)) %>%
      mutate(modulation = x,.before=1)) %>%
  do.call(rbind,.) %>%
  mutate(p.adjust = p.adjust(p.value, method ="bonferroni"), .after = "p.value") %>%
  janitor::clean_names()

ww_tests

## # A tibble: 8 × 6
##   modulation statistic p_value p_adjust parameter method                                       
##   <fct>          <dbl>   <dbl>    <dbl>     <dbl> <chr>                                        
## 1 4Hz           282.     0        0             4 Watson-Wheeler test for homogeneity of angles
## 2 8Hz           274.     0        0             4 Watson-Wheeler test for homogeneity of angles
## 3 16Hz          143.     0        0             4 Watson-Wheeler test for homogeneity of angles
## 4 25Hz           45.4    0        0             4 Watson-Wheeler test for homogeneity of angles
## 5 33Hz           32.7    0        0             4 Watson-Wheeler test for homogeneity of angles
## 6 40Hz           20.6    0        0             4 Watson-Wheeler test for homogeneity of angles
## 7 50Hz            8.78   0.067    0.536         4 Watson-Wheeler test for homogeneity of angles
## 8 80Hz            2.29   0.683    1             4 Watson-Wheeler test for homogeneity of angles

Rao’s test of homogeneity of means, dispersions

Test if means or vector lengths were significantly different between conditions and, separately, across modulation rates:

Compare playback conditions within each modulation rate

How did distributions of call onsets differ between playback conditions at each modulation rate? In their angular means or their polar concentrations?

rao_testsm <- do.call(rbind, lapply(sort(unique(circ_data_df$modulation)), function(x) {
    df <- circ_data_df[circ_data_df$modulation==x,] %>% group_split(condition)
    rao <- circular::rao.test(lapply(df,'[[', 'phase'))
    data.frame(modulation = as.factor(x),
               test = as.factor(c("polar vectors", "dispersions")), 
               statistic = rao$statistic, df = rao$df, p_value = rao$p.value) %>% 
      mutate(across(where(is.numeric), round, 3)) %>% 
      mutate(p_adjust = p.adjust(p_value, method ="bonferroni"), .after = "p_value") 
    }))

# test differences for playback conditions within mod rates
rao_testsm

##    modulation          test statistic df p_value p_adjust
## 1         4Hz polar vectors     3.110  2   0.211    0.422
## 2         4Hz   dispersions   232.370  2   0.000    0.000
## 3         8Hz polar vectors    18.883  2   0.000    0.000
## 4         8Hz   dispersions   199.932  2   0.000    0.000
## 5        16Hz polar vectors     6.895  2   0.032    0.064
## 6        16Hz   dispersions    77.400  2   0.000    0.000
## 7        25Hz polar vectors     3.853  2   0.146    0.292
## 8        25Hz   dispersions    22.588  2   0.000    0.000
## 9        33Hz polar vectors     1.668  2   0.434    0.868
## 10       33Hz   dispersions     7.870  2   0.020    0.040
## 11       40Hz polar vectors     0.897  2   0.639    1.000
## 12       40Hz   dispersions     7.848  2   0.020    0.040
## 13       50Hz polar vectors     0.106  2   0.948    1.000
## 14       50Hz   dispersions     2.521  2   0.284    0.568
## 15       80Hz polar vectors     0.806  2   0.668    1.000
## 16       80Hz   dispersions     0.492  2   0.782    1.000

rao_testsc <- do.call(rbind, lapply(sort(unique(circ_data_df$condition)), function(x) {
    df <- circ_data_df[circ_data_df$condition==x,] %>% group_split(modulation)
    rao <- circular::rao.test(lapply(df,'[[', 'phase'))
    data.frame(condition = as.factor(x),
               test = as.factor(c("polar vectors", "dispersions")), 
               statistic = rao$statistic, df = rao$df, p_value = rao$p.value) %>% 
      mutate(across(where(is.numeric), round, 3)) %>% 
      mutate(p_adjust = p.adjust(p_value, method ="bonferroni"), .after = "p_value")
    }))

# test differences for mod rates within conditions
rao_testsc 

##             condition          test statistic df p_value p_adjust
## 1             silence polar vectors     1.220  7   0.990    1.000
## 2             silence   dispersions     1.016  7   0.995    1.000
## 3 steady-state masker polar vectors    36.305  7   0.000    0.000
## 4 steady-state masker   dispersions   233.993  7   0.000    0.000
## 5       random masker polar vectors    13.431  7   0.062    0.124
## 6       random masker   dispersions   277.880  7   0.000    0.000

Post-hoc

Next, let’s compute post-hoc comparisons for each significant test (by modulation rate) above by comparing each pair of masking conditions within each modulation condition for both rates, again using the Rao test.

do_contrasts <- rao_testsm$modulation[rao_testsm$p_value<0.05]

circ_contrasts <- lapply(do_contrasts, function(x) {
  df <- circ_data_df[circ_data_df$modulation==x,]
  pairs <- t(combn(unique(df$condition), m=2))

  tests <- lapply(seq.int(nrow(pairs)), function(y) {
    condlist <- df[df$condition %in% pairs[y,],] %>% group_split(condition)
    rao <- circular::rao.test(condlist[[1]]$phase, condlist[[2]]$phase)
    data.frame(modulation = as.factor(x),
               condition = as.factor(paste(pairs[y,],collapse = " vs. ")),
               test = as.factor(c("polar vectors", "dispersions")), 
               statistic = rao$statistic, df = rao$df, p_value = rao$p.value) %>% 
      mutate(across(is.numeric, round, 3)) %>% 
      mutate(p_adjust = p.adjust(p_value, method ="bonferroni"), .after = "p_value")
    }) %>% do.call(rbind,.) 
  }
 )  %>% do.call(rbind,.)

circ_contrasts

##    modulation                             condition          test statistic df p_value p_adjust
## 1         4Hz       silence vs. steady-state masker polar vectors     2.985  1   0.084    0.168
## 2         4Hz       silence vs. steady-state masker   dispersions    79.850  1   0.000    0.000
## 3         4Hz             silence vs. random masker polar vectors     3.087  1   0.079    0.158
## 4         4Hz             silence vs. random masker   dispersions   153.051  1   0.000    0.000
## 5         4Hz steady-state masker vs. random masker polar vectors     0.000  1   0.983    1.000
## 6         4Hz steady-state masker vs. random masker   dispersions    15.231  1   0.000    0.000
## 7         8Hz       silence vs. steady-state masker polar vectors     0.046  1   0.830    1.000
## 8         8Hz       silence vs. steady-state masker   dispersions   104.750  1   0.000    0.000
## 9         8Hz             silence vs. random masker polar vectors     0.000  1   0.985    1.000
## 10        8Hz             silence vs. random masker   dispersions    95.520  1   0.000    0.000
## 11        8Hz steady-state masker vs. random masker polar vectors    18.858  1   0.000    0.000
## 12        8Hz steady-state masker vs. random masker   dispersions     0.944  1   0.331    0.662
## 13        8Hz       silence vs. steady-state masker polar vectors     0.046  1   0.830    1.000
## 14        8Hz       silence vs. steady-state masker   dispersions   104.750  1   0.000    0.000
## 15        8Hz             silence vs. random masker polar vectors     0.000  1   0.985    1.000
## 16        8Hz             silence vs. random masker   dispersions    95.520  1   0.000    0.000
## 17        8Hz steady-state masker vs. random masker polar vectors    18.858  1   0.000    0.000
## 18        8Hz steady-state masker vs. random masker   dispersions     0.944  1   0.331    0.662
## 19       16Hz       silence vs. steady-state masker polar vectors     0.848  1   0.357    0.714
## 20       16Hz       silence vs. steady-state masker   dispersions    51.665  1   0.000    0.000
## 21       16Hz             silence vs. random masker polar vectors     1.579  1   0.209    0.418
## 22       16Hz             silence vs. random masker   dispersions    25.990  1   0.000    0.000
## 23       16Hz steady-state masker vs. random masker polar vectors     5.969  1   0.015    0.030
## 24       16Hz steady-state masker vs. random masker   dispersions     0.470  1   0.493    0.986
## 25       16Hz       silence vs. steady-state masker polar vectors     0.848  1   0.357    0.714
## 26       16Hz       silence vs. steady-state masker   dispersions    51.665  1   0.000    0.000
## 27       16Hz             silence vs. random masker polar vectors     1.579  1   0.209    0.418
## 28       16Hz             silence vs. random masker   dispersions    25.990  1   0.000    0.000
## 29       16Hz steady-state masker vs. random masker polar vectors     5.969  1   0.015    0.030
## 30       16Hz steady-state masker vs. random masker   dispersions     0.470  1   0.493    0.986
## 31       25Hz       silence vs. steady-state masker polar vectors     0.044  1   0.834    1.000
## 32       25Hz       silence vs. steady-state masker   dispersions    16.873  1   0.000    0.000
## 33       25Hz             silence vs. random masker polar vectors     0.008  1   0.927    1.000
## 34       25Hz             silence vs. random masker   dispersions     5.805  1   0.016    0.032
## 35       25Hz steady-state masker vs. random masker polar vectors     3.840  1   0.050    0.100
## 36       25Hz steady-state masker vs. random masker   dispersions     0.732  1   0.392    0.784
## 37       33Hz       silence vs. steady-state masker polar vectors     1.651  1   0.199    0.398
## 38       33Hz       silence vs. steady-state masker   dispersions     0.629  1   0.428    0.856
## 39       33Hz             silence vs. random masker polar vectors     0.016  1   0.899    1.000
## 40       33Hz             silence vs. random masker   dispersions     7.305  1   0.007    0.014
## 41       33Hz steady-state masker vs. random masker polar vectors     0.018  1   0.892    1.000
## 42       33Hz steady-state masker vs. random masker   dispersions     6.521  1   0.011    0.022
## 43       40Hz       silence vs. steady-state masker polar vectors     0.849  1   0.357    0.714
## 44       40Hz       silence vs. steady-state masker   dispersions     6.625  1   0.010    0.020
## 45       40Hz             silence vs. random masker polar vectors     0.304  1   0.582    1.000
## 46       40Hz             silence vs. random masker   dispersions     1.231  1   0.267    0.534
## 47       40Hz steady-state masker vs. random masker polar vectors     0.289  1   0.591    1.000
## 48       40Hz steady-state masker vs. random masker   dispersions     0.362  1   0.547    1.000

Rate of vocalizations

Let’s look at the overall number of calls observed in the three experimental conditions and within the two tested modulation rates.

Chi-square test for number of observed calls [not used in paper]

Two post-hoc hypotheses: 1. The number of calls observed in the silent condition would be be greater than that in the masking conditions. - Method: Pearson’s Chi-squared Test for Goodness-of-Fit - H1: Observed calls unevenly distributed across groups.

2. Fewer calls would be observed for the faster masking conditions compared with the slower masking conditions, due to the greater masking effect elicited by more rapid (i.e. “noisier”) amplitude modulations.

• Method: Pearson’s Chi-squared Test for Homogeneity
• H1: Two samples differ in their distribution between groups (levels of masking).

Contingency tables:

n_calls <- data %>% group_by(modulation, condition) %>% summarise(n = n())

(n_conttable <- addmargins(table(data$modulation,data$condition), c(1,2)))

##       
##        silence steady-state masker random masker     Sum
##   4Hz    59635               28938        105990  194563
##   8Hz    62203               65446         51939  179588
##   16Hz   65769               60482         25543  151794
##   25Hz   67403               32743         16600  116746
##   33Hz   78179               34843         12830  125852
##   40Hz   68797               34176         10531  113504
##   50Hz   61956               25425          8487   95868
##   80Hz   75144               30675          5260  111079
##   Sum   539086              312728        237180 1088994

Goodness of fit tests (test of imbalanced observations) for conditions within each modulation rate:

##   4Hz :  X(2, N = 194563) = 46401.95, p < 0.01 (2-tailed) 
##   8Hz :  X(2, N = 179588) = 1661.06, p < 0.01 (2-tailed) 
##   16Hz :  X(2, N = 151794) = 18886.24, p < 0.01 (2-tailed) 
##   25Hz :  X(2, N = 116746) = 34629.51, p < 0.01 (2-tailed) 
##   33Hz :  X(2, N = 125852) = 52705.34, p < 0.01 (2-tailed) 
##   40Hz :  X(2, N = 113504) = 45395.99, p < 0.01 (2-tailed) 
##   50Hz :  X(2, N = 95868) = 46734.51, p < 0.01 (2-tailed) 
##   80Hz :  X(2, N = 111079) = 67584.22, p < 0.01 (2-tailed)

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL
## 
## [[5]]
## NULL
## 
## [[6]]
## NULL
## 
## [[7]]
## NULL
## 
## [[8]]
## NULL

Goodness of fit tests for rates within conditions:

##   silence :  X(7, N = 539086) = 4417.96, p < 0.01 (2-tailed) 
##   steady-state masker :  X(7, N = 312728) = 40810.74, p < 0.01 (2-tailed) 
##   random masker :  X(7, N = 237180) = 276682.97, p < 0.01 (2-tailed)

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL

In the chi-square results below, columns refer to the following data:

• observed: the observed counts,
• expected: the expected counts under the null hypothesis,
• residuals: the Pearson residuals, \((observed - expected) / \sqrt(expected)\),
• stdres: standardized residuals, \((observed - expected) / \sqrt(V)\),

where \(V\) is the residual cell variance.

## Pearson's Chi-squared test for each modulation rate x condition:

##   Mods vs. Conds:  X(14, N = 1088994) = 209611.08, p < 0.01 (2-tailed)

## # A tibble: 24 × 9
##    condition           modulation observed  prop row_prop col_prop expected  resid std_resid
##    <fct>               <fct>         <dbl> <dbl>    <dbl>    <dbl>    <dbl>  <dbl>     <dbl>
##  1 silence             4Hz           59635 0.055    0.111    0.307   96315. -118.     -184. 
##  2 steady-state masker 4Hz           28938 0.027    0.093    0.149   55873. -114.     -149. 
##  3 random masker       4Hz          105990 0.097    0.447    0.545   42375.  309.      386. 
##  4 silence             8Hz           62203 0.057    0.115    0.346   88902.  -89.5    -138. 
##  5 steady-state masker 8Hz           65446 0.06     0.209    0.364   51573.   61.1      79.2
##  6 random masker       8Hz           51939 0.048    0.219    0.289   39114.   64.8      80.2
##  7 silence             16Hz          65769 0.06     0.122    0.433   75143.  -34.2     -51.9
##  8 steady-state masker 16Hz          60482 0.056    0.193    0.398   43591.   80.9     103. 
##  9 random masker       16Hz          25543 0.023    0.108    0.168   33060.  -41.3     -50.4
## 10 silence             25Hz          67403 0.062    0.125    0.577   57793.   40.0      59.5
## # … with 14 more rows

Negative Binomial Regression

First we modelled our counts using a Poisson GLM. However, we observed that the data were highly overdispersed. For example, for the 8 Hz context, overdispersion was:

# Poisson Regression 
calls.glm <- glm(n ~  condition, family = poisson, data = n_calls_reg[n_calls_reg$modulation=="8Hz",])
# model summary
  # summary(calls.glm)

# check for overdispersion
OD.pr <- sum(residuals(calls.glm, type="pearson")^2)/df.residual(calls.glm)
print(OD.pr)

## [1] 6718.167

Therefore, we proceeded with a Negative Binomial GLM, which allows for overdispersion in the count variable.

We ran a model for each modulation rate separately, so that the reference group (Intercept) is meaningful for all model terms for our purposes, i.e. coefficient indicate changes in calling rate between silent baseline and each masking condition.

Our models are given by: \[ln(\widehat{calls_{i}}) = Intercept + \beta_1 Mask(condition_j = 2) + \beta_3 Mask(condition_j = 3)\]

Where \(i\) gives the modulation rate and \(j\) gives the level of the condition manipulation.

Thus, \[\widehat{calls_i} = e^{Intercept + \beta_1 Mask_{condition_j = 2}+...}\] defines the rate of calling per hour (the sampling time for each condition), given an arbitrary combination of modulation rate and masking condition.

We also ran a “big” model that included all interactions, for which coefficients would have a more complex meaning but analyses of variance would provide a global picture.

Run all models

big.model <-  run_negbinom_model(n_calls_reg, 
                                 formula = "n ~ condition * modulation", 
                                 contformula = "~ modulation | condition ")

## [1] "~~~~~~~~~~~~~~  all rates  ~~~~~~~~~~~~~~"
## [1] "Model summary..."
## 
## Call:
## MASS::glm.nb(formula = as.formula(formula), data = dat, init.theta = 0.4410757978, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4987  -1.4098  -0.5586   0.2112   1.8856  
## 
## Coefficients:
##                                             Estimate Std. Error z value             Pr(>|z|)    
## (Intercept)                                  8.00027    0.33671  23.760 < 0.0000000000000002 ***
## conditionsteady-state masker                -0.72309    0.47620  -1.518             0.128903    
## conditionrandom masker                       0.57510    0.47618   1.208             0.227144    
## modulation8Hz                                0.04216    0.47618   0.089             0.929449    
## modulation16Hz                               0.09791    0.47618   0.206             0.837099    
## modulation25Hz                               0.12245    0.47618   0.257             0.797069    
## modulation33Hz                               0.27076    0.47618   0.569             0.569625    
## modulation40Hz                               0.19421    0.48241   0.403             0.687252    
## modulation50Hz                               0.03818    0.47618   0.080             0.936092    
## modulation80Hz                               0.23116    0.47618   0.485             0.627355    
## conditionsteady-state masker:modulation8Hz   0.82520    0.67785   1.217             0.223461    
## conditionrandom masker:modulation8Hz        -0.75544    0.67342  -1.122             0.261954    
## conditionsteady-state masker:modulation16Hz  0.63928    0.67344   0.949             0.342477    
## conditionrandom masker:modulation16Hz       -1.52089    0.67344  -2.258             0.023921 *  
## conditionsteady-state masker:modulation25Hz  0.05238    0.67786   0.077             0.938407    
## conditionrandom masker:modulation25Hz       -1.97639    0.67345  -2.935             0.003339 ** 
## conditionsteady-state masker:modulation33Hz -0.08506    0.67345  -0.126             0.899488    
## conditionrandom masker:modulation33Hz       -2.38232    0.67346  -3.537             0.000404 ***
## conditionsteady-state masker:modulation40Hz  0.02345    0.68225   0.034             0.972580    
## conditionrandom masker:modulation40Hz       -2.50323    0.67789  -3.693             0.000222 ***
## conditionsteady-state masker:modulation50Hz -0.16760    0.67346  -0.249             0.803460    
## conditionrandom masker:modulation50Hz       -2.56299    0.67350  -3.806             0.000142 ***
## conditionsteady-state masker:modulation80Hz -0.17287    0.67345  -0.257             0.797414    
## conditionrandom masker:modulation80Hz       -3.23438    0.67355  -4.802           0.00000157 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.4411) family taken to be 1)
## 
##     Null deviance: 713.90  on 475  degrees of freedom
## Residual deviance: 613.84  on 452  degrees of freedom
## AIC: 7837
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.4411 
##           Std. Err.:  0.0234 
## 
##  2 x log-likelihood:  -7786.9700 
## [1] "Model dispersion..."
## [1] 0.85917
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
## 
## Response: n
##                      LR Chisq Df      Pr(>Chisq)    
## condition              45.676  2 0.0000000001207 ***
## modulation             32.121  7 0.0000385773970 ***
## condition:modulation   42.692 14 0.0000959364219 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Incidence rates..."
##                                                  Estimate         2.5 %       97.5 %
## (Intercept)                                 2981.75000000 1645.30186853 6257.4686715
## conditionsteady-state masker                   0.48525195    0.18758610    1.2552607
## conditionrandom masker                         1.77731198    0.68709741    4.5973654
## modulation8Hz                                  1.04306196    0.40323584    2.6981189
## modulation16Hz                                 1.10285906    0.42635344    2.8527930
## modulation25Hz                                 1.13025908    0.43694631    2.9236672
## modulation33Hz                                 1.31095833    0.50680476    3.3910726
## modulation40Hz                                 1.21435222    0.46526044    3.1943692
## modulation50Hz                                 1.03892010    0.40163459    2.6874054
## modulation80Hz                                 1.26006540    0.48712952    3.2594305
## conditionsteady-state masker:modulation8Hz     2.28234285    0.59826840    8.7411771
## conditionrandom masker:modulation8Hz           0.46980603    0.12401175    1.7798129
## conditionsteady-state masker:modulation16Hz    1.89512393    0.50022870    7.1797054
## conditionrandom masker:modulation16Hz          0.21851789    0.05767924    0.8278554
## conditionsteady-state masker:modulation25Hz    1.05377611    0.27621988    4.0359513
## conditionrandom masker:modulation25Hz          0.13856872    0.03657505    0.5249832
## conditionsteady-state masker:modulation33Hz    0.91845555    0.24242820    3.4796307
## conditionrandom masker:modulation33Hz          0.09233639    0.02437154    0.3498346
## conditionsteady-state masker:modulation40Hz    1.02372767    0.26545814    3.9479609
## conditionrandom masker:modulation40Hz          0.08182011    0.02136184    0.3121601
## conditionsteady-state masker:modulation50Hz    0.84568826    0.22321671    3.2040103
## conditionrandom masker:modulation50Hz          0.07707387    0.02034186    0.2920274
## conditionsteady-state masker:modulation80Hz    0.84124592    0.22204717    3.1871367
## conditionrandom masker:modulation80Hz          0.03938472    0.01039366    0.1492406
## [1] "Estimated marginal mean contrasts..."
## $emmeans
## condition = silence:
##  modulation response     SE  df asymp.LCL asymp.UCL
##  4Hz            2982 1004.0 Inf      1541      5769
##  8Hz            3110 1047.2 Inf      1608      6017
##  16Hz           3288 1107.3 Inf      1700      6362
##  25Hz           3370 1134.8 Inf      1742      6520
##  33Hz           3909 1316.2 Inf      2020      7563
##  40Hz           3621 1250.9 Inf      1840      7126
##  50Hz           3098 1043.1 Inf      1601      5993
##  80Hz           3757 1265.1 Inf      1942      7269
## 
## condition = steady-state masker:
##  modulation response     SE  df asymp.LCL asymp.UCL
##  4Hz            1447  487.2 Inf       748      2799
##  8Hz            3445 1189.9 Inf      1750      6779
##  16Hz           3024 1018.3 Inf      1563      5851
##  25Hz           1723  595.4 Inf       876      3392
##  33Hz           1742  586.6 Inf       900      3371
##  40Hz           1799  621.4 Inf       914      3540
##  50Hz           1271  428.1 Inf       657      2460
##  80Hz           1534  516.5 Inf       793      2967
## 
## condition = random masker:
##  modulation response     SE  df asymp.LCL asymp.UCL
##  4Hz            5300 1784.4 Inf      2739     10253
##  8Hz            2597  874.4 Inf      1342      5024
##  16Hz           1277  430.1 Inf       660      2471
##  25Hz            830  279.5 Inf       429      1606
##  33Hz            642  216.1 Inf       332      1241
##  40Hz            527  177.4 Inf       272      1019
##  50Hz            424  142.9 Inf       219       821
##  80Hz            263   88.6 Inf       136       509
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
## condition = silence:
##  contrast     ratio    SE  df null z.ratio p.value
##  4Hz / 8Hz    0.959 0.457 Inf    1  -0.089  1.0000
##  4Hz / 16Hz   0.907 0.432 Inf    1  -0.206  1.0000
##  4Hz / 25Hz   0.885 0.421 Inf    1  -0.257  1.0000
##  4Hz / 33Hz   0.763 0.363 Inf    1  -0.569  1.0000
##  4Hz / 40Hz   0.823 0.397 Inf    1  -0.403  1.0000
##  4Hz / 50Hz   0.963 0.458 Inf    1  -0.080  1.0000
##  4Hz / 80Hz   0.794 0.378 Inf    1  -0.485  1.0000
##  8Hz / 16Hz   0.946 0.450 Inf    1  -0.117  1.0000
##  8Hz / 25Hz   0.923 0.439 Inf    1  -0.169  1.0000
##  8Hz / 33Hz   0.796 0.379 Inf    1  -0.480  1.0000
##  8Hz / 40Hz   0.859 0.414 Inf    1  -0.315  1.0000
##  8Hz / 50Hz   1.004 0.478 Inf    1   0.008  1.0000
##  8Hz / 80Hz   0.828 0.394 Inf    1  -0.397  1.0000
##  16Hz / 25Hz  0.976 0.465 Inf    1  -0.052  1.0000
##  16Hz / 33Hz  0.841 0.401 Inf    1  -0.363  1.0000
##  16Hz / 40Hz  0.908 0.438 Inf    1  -0.200  1.0000
##  16Hz / 50Hz  1.062 0.505 Inf    1   0.125  1.0000
##  16Hz / 80Hz  0.875 0.417 Inf    1  -0.280  1.0000
##  25Hz / 33Hz  0.862 0.411 Inf    1  -0.311  1.0000
##  25Hz / 40Hz  0.931 0.449 Inf    1  -0.149  1.0000
##  25Hz / 50Hz  1.088 0.518 Inf    1   0.177  1.0000
##  25Hz / 80Hz  0.897 0.427 Inf    1  -0.228  1.0000
##  33Hz / 40Hz  1.080 0.521 Inf    1   0.159  1.0000
##  33Hz / 50Hz  1.262 0.601 Inf    1   0.488  1.0000
##  33Hz / 80Hz  1.040 0.495 Inf    1   0.083  1.0000
##  40Hz / 50Hz  1.169 0.564 Inf    1   0.323  1.0000
##  40Hz / 80Hz  0.964 0.465 Inf    1  -0.077  1.0000
##  50Hz / 80Hz  0.824 0.393 Inf    1  -0.405  1.0000
## 
## condition = steady-state masker:
##  contrast     ratio    SE  df null z.ratio p.value
##  4Hz / 8Hz    0.420 0.203 Inf    1  -1.798  1.0000
##  4Hz / 16Hz   0.478 0.228 Inf    1  -1.548  1.0000
##  4Hz / 25Hz   0.840 0.405 Inf    1  -0.362  1.0000
##  4Hz / 33Hz   0.831 0.396 Inf    1  -0.390  1.0000
##  4Hz / 40Hz   0.804 0.388 Inf    1  -0.451  1.0000
##  4Hz / 50Hz   1.138 0.542 Inf    1   0.272  1.0000
##  4Hz / 80Hz   0.943 0.449 Inf    1  -0.122  1.0000
##  8Hz / 16Hz   1.139 0.549 Inf    1   0.270  1.0000
##  8Hz / 25Hz   1.999 0.977 Inf    1   1.417  1.0000
##  8Hz / 33Hz   1.977 0.954 Inf    1   1.413  1.0000
##  8Hz / 40Hz   1.915 0.936 Inf    1   1.330  1.0000
##  8Hz / 50Hz   2.710 1.307 Inf    1   2.066  1.0000
##  8Hz / 80Hz   2.246 1.083 Inf    1   1.677  1.0000
##  16Hz / 25Hz  1.755 0.847 Inf    1   1.166  1.0000
##  16Hz / 33Hz  1.736 0.827 Inf    1   1.158  1.0000
##  16Hz / 40Hz  1.681 0.811 Inf    1   1.077  1.0000
##  16Hz / 50Hz  2.379 1.133 Inf    1   1.820  1.0000
##  16Hz / 80Hz  1.972 0.939 Inf    1   1.426  1.0000
##  25Hz / 33Hz  0.989 0.477 Inf    1  -0.023  1.0000
##  25Hz / 40Hz  0.958 0.468 Inf    1  -0.088  1.0000
##  25Hz / 50Hz  1.356 0.654 Inf    1   0.631  1.0000
##  25Hz / 80Hz  1.124 0.542 Inf    1   0.242  1.0000
##  33Hz / 40Hz  0.969 0.467 Inf    1  -0.066  1.0000
##  33Hz / 50Hz  1.370 0.653 Inf    1   0.662  1.0000
##  33Hz / 80Hz  1.136 0.541 Inf    1   0.268  1.0000
##  40Hz / 50Hz  1.415 0.683 Inf    1   0.719  1.0000
##  40Hz / 80Hz  1.173 0.566 Inf    1   0.330  1.0000
##  50Hz / 80Hz  0.829 0.395 Inf    1  -0.394  1.0000
## 
## condition = random masker:
##  contrast     ratio    SE  df null z.ratio p.value
##  4Hz / 8Hz    2.041 0.972 Inf    1   1.498  1.0000
##  4Hz / 16Hz   4.149 1.976 Inf    1   2.988  0.0786
##  4Hz / 25Hz   6.385 3.041 Inf    1   3.893  0.0028
##  4Hz / 33Hz   8.261 3.934 Inf    1   4.434  0.0003
##  4Hz / 40Hz  10.065 4.793 Inf    1   4.848  <.0001
##  4Hz / 50Hz  12.489 5.948 Inf    1   5.301  <.0001
##  4Hz / 80Hz  20.150 9.599 Inf    1   6.305  <.0001
##  8Hz / 16Hz   2.033 0.968 Inf    1   1.490  1.0000
##  8Hz / 25Hz   3.129 1.490 Inf    1   2.395  0.4651
##  8Hz / 33Hz   4.048 1.928 Inf    1   2.936  0.0931
##  8Hz / 40Hz   4.932 2.349 Inf    1   3.351  0.0226
##  8Hz / 50Hz   6.120 2.915 Inf    1   3.803  0.0040
##  8Hz / 80Hz   9.874 4.704 Inf    1   4.807  <.0001
##  16Hz / 25Hz  1.539 0.733 Inf    1   0.905  1.0000
##  16Hz / 33Hz  1.991 0.948 Inf    1   1.446  1.0000
##  16Hz / 40Hz  2.426 1.155 Inf    1   1.860  1.0000
##  16Hz / 50Hz  3.010 1.434 Inf    1   2.313  0.5799
##  16Hz / 80Hz  4.856 2.313 Inf    1   3.317  0.0255
##  25Hz / 33Hz  1.294 0.616 Inf    1   0.541  1.0000
##  25Hz / 40Hz  1.576 0.751 Inf    1   0.955  1.0000
##  25Hz / 50Hz  1.956 0.932 Inf    1   1.408  1.0000
##  25Hz / 80Hz  3.156 1.504 Inf    1   2.412  0.4438
##  33Hz / 40Hz  1.218 0.580 Inf    1   0.415  1.0000
##  33Hz / 50Hz  1.512 0.720 Inf    1   0.868  1.0000
##  33Hz / 80Hz  2.439 1.162 Inf    1   1.872  1.0000
##  40Hz / 50Hz  1.241 0.591 Inf    1   0.453  1.0000
##  40Hz / 80Hz  2.002 0.954 Inf    1   1.457  1.0000
##  50Hz / 80Hz  1.613 0.769 Inf    1   1.004  1.0000
## 
## P value adjustment: bonferroni method for 28 tests 
## Tests are performed on the log scale

# overdispersion
print(big.model$dispersion)

## [1] 0.85917

all.models <- lapply(sort(unique(n_calls_reg$modulation)), function(x) 
                      md <- run_negbinom_model(n_calls_reg %>% filter(modulation==x)) )

## [1] "~~~~~~~~~~~~~~  4Hz  ~~~~~~~~~~~~~~"
## [1] "Model summary..."
## 
## Call:
## MASS::glm.nb(formula = n ~ condition, data = dat, init.theta = 0.4723339773, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3464  -1.3494  -0.7116   0.1990   1.5790  
## 
## Coefficients:
##                              Estimate Std. Error z value            Pr(>|z|)    
## (Intercept)                    8.0003     0.3254  24.587 <0.0000000000000002 ***
## conditionsteady-state masker  -0.7231     0.4602  -1.571               0.116    
## conditionrandom masker         0.5751     0.4602   1.250               0.211    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.4723) family taken to be 1)
## 
##     Null deviance: 84.212  on 59  degrees of freedom
## Residual deviance: 76.661  on 57  degrees of freedom
## AIC: 1051.1
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.4723 
##           Std. Err.:  0.0710 
## 
##  2 x log-likelihood:  -1043.1210 
## [1] "Model dispersion..."
## [1] 0.9538492
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
## 
## Response: n
##           LR Chisq Df Pr(>Chisq)  
## condition   7.5502  2    0.02294 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Incidence rates..."
##                                  Estimate        2.5 %      97.5 %
## (Intercept)                  2981.7500000 1675.1885512 6085.738182
## conditionsteady-state masker    0.4852519    0.1938973    1.214403
## conditionrandom masker          1.7773120    0.7102157    4.447716
## [1] "Estimated marginal mean contrasts..."
## $emmeans
##  condition           response   SE  df asymp.LCL asymp.UCL
##  silence                 2982  970 Inf      1576      5642
##  steady-state masker     1447  471 Inf       765      2738
##  random masker           5300 1724 Inf      2801     10027
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                              ratio    SE  df null z.ratio p.value
##  silence / (steady-state masker)       2.061 0.948 Inf    1   1.571  0.3483
##  silence / random masker               0.563 0.259 Inf    1  -1.250  0.6341
##  (steady-state masker) / random masker 0.273 0.126 Inf    1  -2.821  0.0144
## 
## P value adjustment: bonferroni method for 3 tests 
## Tests are performed on the log scale 
## 
## [1] "~~~~~~~~~~~~~~  8Hz  ~~~~~~~~~~~~~~"
## [1] "Model summary..."
## 
## Call:
## MASS::glm.nb(formula = n ~ condition, data = dat, init.theta = 0.3868466812, 
##     link = log)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.32263  -1.40195  -0.71589   0.06566   1.42301  
## 
## Coefficients:
##                              Estimate Std. Error z value            Pr(>|z|)    
## (Intercept)                    8.0424     0.3595  22.369 <0.0000000000000002 ***
## conditionsteady-state masker   0.1021     0.5151   0.198               0.843    
## conditionrandom masker        -0.1803     0.5085  -0.355               0.723    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.3868) family taken to be 1)
## 
##     Null deviance: 77.915  on 58  degrees of freedom
## Residual deviance: 77.609  on 56  degrees of freedom
## AIC: 1017.8
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.3868 
##           Std. Err.:  0.0577 
## 
##  2 x log-likelihood:  -1009.7670 
## [1] "Model dispersion..."
## [1] 0.8545476
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
## 
## Response: n
##           LR Chisq Df Pr(>Chisq)
## condition   0.3059  2     0.8582
## [1] "Incidence rates..."
##                                  Estimate        2.5 %      97.5 %
## (Intercept)                  3110.1500000 1655.2222327 6909.314899
## conditionsteady-state masker    1.1075113    0.3966732    3.119991
## conditionrandom masker          0.8349919    0.3018830    2.309542
## [1] "Estimated marginal mean contrasts..."
## $emmeans
##  condition           response   SE  df asymp.LCL asymp.UCL
##  silence                 3110 1118 Inf      1537      6292
##  steady-state masker     3445 1271 Inf      1672      7098
##  random masker           2597  934 Inf      1284      5254
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                              ratio    SE  df null z.ratio p.value
##  silence / (steady-state masker)       0.903 0.465 Inf    1  -0.198  1.0000
##  silence / random masker               1.198 0.609 Inf    1   0.355  1.0000
##  (steady-state masker) / random masker 1.326 0.683 Inf    1   0.548  1.0000
## 
## P value adjustment: bonferroni method for 3 tests 
## Tests are performed on the log scale 
## 
## [1] "~~~~~~~~~~~~~~  16Hz  ~~~~~~~~~~~~~~"
## [1] "Model summary..."
## 
## Call:
## MASS::glm.nb(formula = n ~ condition, data = dat, init.theta = 0.5221485665, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2929  -1.2910  -0.4686   0.3407   1.5296  
## 
## Coefficients:
##                              Estimate Std. Error z value            Pr(>|z|)    
## (Intercept)                    8.0982     0.3095  26.168 <0.0000000000000002 ***
## conditionsteady-state masker  -0.0838     0.4377  -0.191              0.8482    
## conditionrandom masker        -0.9458     0.4377  -2.161              0.0307 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.5221) family taken to be 1)
## 
##     Null deviance: 80.629  on 59  degrees of freedom
## Residual deviance: 75.558  on 57  degrees of freedom
## AIC: 1036.4
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.5221 
##           Std. Err.:  0.0792 
## 
##  2 x log-likelihood:  -1028.4460 
## [1] "Model dispersion..."
## [1] 0.8309316
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
## 
## Response: n
##           LR Chisq Df Pr(>Chisq)  
## condition   5.0715  2     0.0792 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Incidence rates..."
##                                  Estimate        2.5 %       97.5 %
## (Intercept)                  3288.4500000 1895.4786790 6456.3843402
## conditionsteady-state masker    0.9196126    0.3848699    2.1973327
## conditionrandom masker          0.3883745    0.1625319    0.9280312
## [1] "Estimated marginal mean contrasts..."
## $emmeans
##  condition           response   SE  df asymp.LCL asymp.UCL
##  silence                 3288 1018 Inf      1793      6031
##  steady-state masker     3024  936 Inf      1649      5547
##  random masker           1277  395 Inf       696      2343
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                              ratio    SE  df null z.ratio p.value
##  silence / (steady-state masker)        1.09 0.476 Inf    1   0.191  1.0000
##  silence / random masker                2.57 1.127 Inf    1   2.161  0.0921
##  (steady-state masker) / random masker  2.37 1.036 Inf    1   1.969  0.1467
## 
## P value adjustment: bonferroni method for 3 tests 
## Tests are performed on the log scale 
## 
## [1] "~~~~~~~~~~~~~~  25Hz  ~~~~~~~~~~~~~~"
## [1] "Model summary..."
## 
## Call:
## MASS::glm.nb(formula = n ~ condition, data = dat, init.theta = 0.4039458612, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2513  -1.3790  -0.7238   0.1119   1.3833  
## 
## Coefficients:
##                              Estimate Std. Error z value             Pr(>|z|)    
## (Intercept)                    8.1227     0.3518  23.086 < 0.0000000000000002 ***
## conditionsteady-state masker  -0.6707     0.5041  -1.331              0.18335    
## conditionrandom masker        -1.4013     0.4976  -2.816              0.00486 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.4039) family taken to be 1)
## 
##     Null deviance: 84.625  on 58  degrees of freedom
## Residual deviance: 77.072  on 56  degrees of freedom
## AIC: 954.96
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.4039 
##           Std. Err.:  0.0604 
## 
##  2 x log-likelihood:  -946.9580 
## [1] "Model dispersion..."
## [1] 0.8618375
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
## 
## Response: n
##           LR Chisq Df Pr(>Chisq)  
## condition   7.5536  2     0.0229 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Incidence rates..."
##                                  Estimate         2.5 %       97.5 %
## (Intercept)                  3370.1500000 1815.42972059 7344.0187496
## conditionsteady-state masker    0.5113469    0.18736028    1.4075862
## conditionrandom masker          0.2462798    0.09107138    0.6660024
## [1] "Estimated marginal mean contrasts..."
## $emmeans
##  condition           response   SE  df asymp.LCL asymp.UCL
##  silence                 3370 1186 Inf      1691      6716
##  steady-state masker     1723  622 Inf       849      3497
##  random masker            830  292 Inf       416      1654
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                              ratio    SE  df null z.ratio p.value
##  silence / (steady-state masker)        1.96 0.986 Inf    1   1.331  0.5501
##  silence / random masker                4.06 2.021 Inf    1   2.816  0.0146
##  (steady-state masker) / random masker  2.08 1.047 Inf    1   1.449  0.4419
## 
## P value adjustment: bonferroni method for 3 tests 
## Tests are performed on the log scale 
## 
## [1] "~~~~~~~~~~~~~~  33Hz  ~~~~~~~~~~~~~~"
## [1] "Model summary..."
## 
## Call:
## MASS::glm.nb(formula = n ~ condition, data = dat, init.theta = 0.4259724824, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4566  -1.3961  -0.5122   0.2072   1.8530  
## 
## Coefficients:
##                              Estimate Std. Error z value             Pr(>|z|)    
## (Intercept)                    8.2710     0.3426  24.140 < 0.0000000000000002 ***
## conditionsteady-state masker  -0.8081     0.4846  -1.668             0.095356 .  
## conditionrandom masker        -1.8072     0.4846  -3.729             0.000192 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.426) family taken to be 1)
## 
##     Null deviance: 90.429  on 59  degrees of freedom
## Residual deviance: 77.690  on 57  degrees of freedom
## AIC: 973.91
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.4260 
##           Std. Err.:  0.0635 
## 
##  2 x log-likelihood:  -965.9070 
## [1] "Model dispersion..."
## [1] 0.8952745
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
## 
## Response: n
##           LR Chisq Df Pr(>Chisq)   
## condition   12.739  2   0.001713 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Incidence rates..."
##                                  Estimate         2.5 %       97.5 %
## (Intercept)                  3908.9500000 2136.57666836 8324.8275436
## conditionsteady-state masker    0.4456823    0.16933647    1.1730063
## conditionrandom masker          0.1641106    0.06234774    0.4319688
## [1] "Estimated marginal mean contrasts..."
## $emmeans
##  condition           response   SE  df asymp.LCL asymp.UCL
##  silence                 3909 1339 Inf      1997      7651
##  steady-state masker     1742  597 Inf       890      3410
##  random masker            642  220 Inf       328      1256
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                              ratio   SE  df null z.ratio p.value
##  silence / (steady-state masker)        2.24 1.09 Inf    1   1.668  0.2861
##  silence / random masker                6.09 2.95 Inf    1   3.729  0.0006
##  (steady-state masker) / random masker  2.72 1.32 Inf    1   2.062  0.1178
## 
## P value adjustment: bonferroni method for 3 tests 
## Tests are performed on the log scale 
## 
## [1] "~~~~~~~~~~~~~~  40Hz  ~~~~~~~~~~~~~~"
## [1] "Model summary..."
## 
## Call:
## MASS::glm.nb(formula = n ~ condition, data = dat, init.theta = 0.4689172885, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0805  -1.4607  -0.5346   0.3031   1.4668  
## 
## Coefficients:
##                              Estimate Std. Error z value             Pr(>|z|)    
## (Intercept)                    8.1945     0.3350  24.458 < 0.0000000000000002 ***
## conditionsteady-state masker  -0.6996     0.4738  -1.477                 0.14    
## conditionrandom masker        -1.9281     0.4680  -4.120            0.0000378 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.4689) family taken to be 1)
## 
##     Null deviance: 89.353  on 57  degrees of freedom
## Residual deviance: 74.197  on 55  degrees of freedom
## AIC: 940.19
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.4689 
##           Std. Err.:  0.0716 
## 
##  2 x log-likelihood:  -932.1910 
## [1] "Model dispersion..."
## [1] 0.841219
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
## 
## Response: n
##           LR Chisq Df Pr(>Chisq)    
## condition   15.156  2  0.0005115 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Incidence rates..."
##                                  Estimate         2.5 %       97.5 %
## (Intercept)                  3620.8947368 2003.10860148 7568.0222469
## conditionsteady-state masker    0.4967658    0.19297987    1.2787671
## conditionrandom masker          0.1454199    0.05697374    0.3684498
## [1] "Estimated marginal mean contrasts..."
## $emmeans
##  condition           response   SE  df asymp.LCL asymp.UCL
##  silence                 3621 1213 Inf      1878      6982
##  steady-state masker     1799  603 Inf       933      3469
##  random masker            527  172 Inf       278       999
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                              ratio    SE  df null z.ratio p.value
##  silence / (steady-state masker)        2.01 0.954 Inf    1   1.477  0.4194
##  silence / random masker                6.88 3.218 Inf    1   4.120  0.0001
##  (steady-state masker) / random masker  3.42 1.599 Inf    1   2.625  0.0260
## 
## P value adjustment: bonferroni method for 3 tests 
## Tests are performed on the log scale 
## 
## [1] "~~~~~~~~~~~~~~  50Hz  ~~~~~~~~~~~~~~"
## [1] "Model summary..."
## 
## Call:
## MASS::glm.nb(formula = n ~ condition, data = dat, init.theta = 0.4025397372, 
##     link = log)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.17513  -1.47676  -0.69721   0.04592   1.42237  
## 
## Coefficients:
##                              Estimate Std. Error z value             Pr(>|z|)    
## (Intercept)                    8.0384     0.3525  22.807 < 0.0000000000000002 ***
## conditionsteady-state masker  -0.8907     0.4985  -1.787                0.074 .  
## conditionrandom masker        -1.9879     0.4986  -3.987            0.0000668 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.4025) family taken to be 1)
## 
##     Null deviance: 92.781  on 59  degrees of freedom
## Residual deviance: 78.417  on 57  degrees of freedom
## AIC: 928.2
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.4025 
##           Std. Err.:  0.0597 
## 
##  2 x log-likelihood:  -920.2040 
## [1] "Model dispersion..."
## [1] 0.8204445
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
## 
## Response: n
##           LR Chisq Df Pr(>Chisq)    
## condition   14.364  2  0.0007603 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Incidence rates..."
##                                  Estimate         2.5 %       97.5 %
## (Intercept)                  3097.8000000 1667.09359936 6760.9100740
## conditionsteady-state masker    0.4103719    0.15148246    1.1117133
## conditionrandom masker          0.1369843    0.05055829    0.3711498
## [1] "Estimated marginal mean contrasts..."
## $emmeans
##  condition           response   SE  df asymp.LCL asymp.UCL
##  silence                 3098 1092 Inf      1553      6181
##  steady-state masker     1271  448 Inf       637      2537
##  random masker            424  150 Inf       213       847
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                              ratio   SE  df null z.ratio p.value
##  silence / (steady-state masker)        2.44 1.21 Inf    1   1.787  0.2219
##  silence / random masker                7.30 3.64 Inf    1   3.987  0.0002
##  (steady-state masker) / random masker  3.00 1.49 Inf    1   2.201  0.0833
## 
## P value adjustment: bonferroni method for 3 tests 
## Tests are performed on the log scale 
## 
## [1] "~~~~~~~~~~~~~~  80Hz  ~~~~~~~~~~~~~~"
## [1] "Model summary..."
## 
## Call:
## MASS::glm.nb(formula = n ~ condition, data = dat, init.theta = 0.4827242071, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.0092  -1.3175  -0.5869   0.2676   1.4456  
## 
## Coefficients:
##                              Estimate Std. Error z value             Pr(>|z|)    
## (Intercept)                    8.2314     0.3219  25.575 < 0.0000000000000002 ***
## conditionsteady-state masker  -0.8960     0.4552  -1.968                0.049 *  
## conditionrandom masker        -2.6593     0.4554  -5.840        0.00000000523 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.4827) family taken to be 1)
## 
##     Null deviance: 103.877  on 59  degrees of freedom
## Residual deviance:  76.242  on 57  degrees of freedom
## AIC: 944.88
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.4827 
##           Std. Err.:  0.0728 
## 
##  2 x log-likelihood:  -936.8780 
## [1] "Model dispersion..."
## [1] 0.8198301
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
## 
## Response: n
##           LR Chisq Df   Pr(>Chisq)    
## condition   27.634  2 0.0000009984 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Incidence rates..."
##                                   Estimate        2.5 %       97.5 %
## (Intercept)                  3757.20000000 2122.8400904 7602.5475769
## conditionsteady-state masker    0.40821622    0.1647983    1.0111782
## conditionrandom masker          0.06999894    0.0282497    0.1734479
## [1] "Estimated marginal mean contrasts..."
## $emmeans
##  condition           response     SE  df asymp.LCL asymp.UCL
##  silence                 3757 1209.3 Inf      1999      7060
##  steady-state masker     1534  493.7 Inf       816      2882
##  random masker            263   84.7 Inf       140       494
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast                              ratio   SE  df null z.ratio p.value
##  silence / (steady-state masker)        2.45 1.12 Inf    1   1.968  0.1471
##  silence / random masker               14.29 6.51 Inf    1   5.840  <.0001
##  (steady-state masker) / random masker  5.83 2.66 Inf    1   3.872  0.0003
## 
## P value adjustment: bonferroni method for 3 tests 
## Tests are performed on the log scale

names(all.models) <- sort(unique(n_calls_reg$modulation))

# overdispersion
print(do.call(rbind, lapply(all.models, '[[', 'dispersion')) )

##           [,1]
## 4Hz  0.9538492
## 8Hz  0.8545476
## 16Hz 0.8309316
## 25Hz 0.8618375
## 33Hz 0.8952745
## 40Hz 0.8412190
## 50Hz 0.8204445
## 80Hz 0.8198301

Anova

Note that these are 1-way ANOVAs from different models put together in one table, not a single ANOVA.

# pull out anova tables
(anova <- do.call(rbind, lapply(all.models, '[[', 'deviance')) %>% 
   mutate(across(is.numeric, round, 2)))

##      LR Chisq Df Pr(>Chisq)
## 4Hz      7.55  2       0.02
## 8Hz      0.31  2       0.86
## 16Hz     5.07  2       0.08
## 25Hz     7.55  2       0.02
## 33Hz    12.74  2       0.00
## 40Hz    15.16  2       0.00
## 50Hz    14.36  2       0.00
## 80Hz    27.63  2       0.00

Below, the analysis of variance for the large model with all features:

big.model$deviance

##                      LR Chisq Df Pr(>Chisq)
## condition              45.676  2          0
## modulation             32.121  7          0
## condition:modulation   42.692 14          0

Incidence Rate Ratios (IRR)

The IRRs give us the change in the per hour rate [of calling] from baseline to each condition.

# IRR
IRRs <- do.call(rbind, lapply(all.models, '[[', 'IRR')) %>%
   mutate(across(is.numeric, round, 2)) %>%
   mutate(modulation = str_split(row.names(.),"[.]",simplify = T)[,1], .before=1) %>%
   mutate(coef = str_remove(str_split(row.names(.),"[.]",simplify = T)[,2],"condition"), .after=1) %>%
   remove_rownames()
IRRs

##    modulation                coef Estimate   2.5 %  97.5 %
## 1         4Hz         (Intercept)  2981.75 1675.19 6085.74
## 2         4Hz steady-state masker     0.49    0.19    1.21
## 3         4Hz       random masker     1.78    0.71    4.45
## 4         8Hz         (Intercept)  3110.15 1655.22 6909.31
## 5         8Hz steady-state masker     1.11    0.40    3.12
## 6         8Hz       random masker     0.83    0.30    2.31
## 7        16Hz         (Intercept)  3288.45 1895.48 6456.38
## 8        16Hz steady-state masker     0.92    0.38    2.20
## 9        16Hz       random masker     0.39    0.16    0.93
## 10       25Hz         (Intercept)  3370.15 1815.43 7344.02
## 11       25Hz steady-state masker     0.51    0.19    1.41
## 12       25Hz       random masker     0.25    0.09    0.67
## 13       33Hz         (Intercept)  3908.95 2136.58 8324.83
## 14       33Hz steady-state masker     0.45    0.17    1.17
## 15       33Hz       random masker     0.16    0.06    0.43
## 16       40Hz         (Intercept)  3620.89 2003.11 7568.02
## 17       40Hz steady-state masker     0.50    0.19    1.28
## 18       40Hz       random masker     0.15    0.06    0.37
## 19       50Hz         (Intercept)  3097.80 1667.09 6760.91
## 20       50Hz steady-state masker     0.41    0.15    1.11
## 21       50Hz       random masker     0.14    0.05    0.37
## 22       80Hz         (Intercept)  3757.20 2122.84 7602.55
## 23       80Hz steady-state masker     0.41    0.16    1.01
## 24       80Hz       random masker     0.07    0.03    0.17

Marginal mean contrasts

Below are the post-hoc contrasts from estimated marginal means:

(contrasts <- do.call(rbind, lapply(all.models, '[[', 'contrasts')) %>%
   mutate(across(is.numeric, round, 2)) %>%
   mutate(modulation = str_split(row.names(.),"[.]",simplify = T)[,1], .before=1) %>%
   remove_rownames() %>%
   dplyr::select(-c("df", "null")) %>% 
   dplyr::rename(IRR = ratio))

##    modulation                     masking_condition   IRR   se asymp_lcl asymp_ucl z_ratio p_value
## 1         4Hz       silence / (steady-state masker)  2.06 0.95      0.69      6.20    1.57    0.35
## 2         4Hz               silence / random masker  0.56 0.26      0.19      1.69   -1.25    0.63
## 3         4Hz (steady-state masker) / random masker  0.27 0.13      0.09      0.82   -2.82    0.01
## 4         8Hz       silence / (steady-state masker)  0.90 0.47      0.26      3.10   -0.20    1.00
## 5         8Hz               silence / random masker  1.20 0.61      0.36      4.04    0.36    1.00
## 6         8Hz (steady-state masker) / random masker  1.33 0.68      0.39      4.55    0.55    1.00
## 7        16Hz       silence / (steady-state masker)  1.09 0.48      0.38      3.10    0.19    1.00
## 8        16Hz               silence / random masker  2.58 1.13      0.90      7.34    2.16    0.09
## 9        16Hz (steady-state masker) / random masker  2.37 1.04      0.83      6.75    1.97    0.15
## 10       25Hz       silence / (steady-state masker)  1.96 0.99      0.58      6.54    1.33    0.55
## 11       25Hz               silence / random masker  4.06 2.02      1.23     13.36    2.82    0.01
## 12       25Hz (steady-state masker) / random masker  2.08 1.05      0.62      6.94    1.45    0.44
## 13       33Hz       silence / (steady-state masker)  2.24 1.09      0.70      7.16    1.67    0.29
## 14       33Hz               silence / random masker  6.09 2.95      1.91     19.44    3.73    0.00
## 15       33Hz (steady-state masker) / random masker  2.72 1.32      0.85      8.66    2.06    0.12
## 16       40Hz       silence / (steady-state masker)  2.01 0.95      0.65      6.26    1.48    0.42
## 17       40Hz               silence / random masker  6.88 3.22      2.24     21.08    4.12    0.00
## 18       40Hz (steady-state masker) / random masker  3.42 1.60      1.11     10.47    2.62    0.03
## 19       50Hz       silence / (steady-state masker)  2.44 1.22      0.74      8.04    1.79    0.22
## 20       50Hz               silence / random masker  7.30 3.64      2.21     24.08    3.99    0.00
## 21       50Hz (steady-state masker) / random masker  3.00 1.49      0.91      9.88    2.20    0.08
## 22       80Hz       silence / (steady-state masker)  2.45 1.11      0.82      7.28    1.97    0.15
## 23       80Hz               silence / random masker 14.29 6.50      4.80     42.50    5.84    0.00
## 24       80Hz (steady-state masker) / random masker  5.83 2.66      1.96     17.35    3.87    0.00

Model predictions

Now let’s plot the model predictions to ensure it’s a good fit:

preds <- do.call(rbind, lapply(all.models, '[[', 'predictions')) %>%
   mutate(across(is.numeric, round, 2)) %>%
   mutate(modulation = factor(str_split(row.names(.),"[.]",simplify = T)[,1], 
          levels = sort(unique(data$modulation))),.before=1) %>%
   remove_rownames()

n_pred <- p.n_pred(preds) %>% recolor()

n_pred

ggsave(file.path(figures_path,paste0(exp, '_fig1e.png') ), n_pred, width = 5, height = 4)

Changes in call rate between modulation rate contexts

Now let’s see if there were changes in the overall rate of calling between different modulation rate contexts for the same playback condition.

all.models2 <- lapply(sort(unique(n_calls_reg$condition)), function(x) 
                      md <- run_negbinom_model(n_calls_reg %>% filter(condition==x), 
                                               formula = "n ~  modulation", 
                                               contformula = "~ modulation")  )

## [1] "~~~~~~~~~~~~~~  all rates  ~~~~~~~~~~~~~~"
## [1] "Model summary..."
## 
## Call:
## MASS::glm.nb(formula = as.formula(formula), data = dat, init.theta = 0.4097400923, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4104  -1.3337  -0.6802   0.2266   1.4645  
## 
## Coefficients:
##                Estimate Std. Error z value            Pr(>|z|)    
## (Intercept)     8.00027    0.34935  22.900 <0.0000000000000002 ***
## modulation8Hz   0.04216    0.49405   0.085               0.932    
## modulation16Hz  0.09791    0.49405   0.198               0.843    
## modulation25Hz  0.12245    0.49405   0.248               0.804    
## modulation33Hz  0.27076    0.49405   0.548               0.584    
## modulation40Hz  0.19421    0.50051   0.388               0.698    
## modulation50Hz  0.03818    0.49405   0.077               0.938    
## modulation80Hz  0.23116    0.49405   0.468               0.640    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.4097) family taken to be 1)
## 
##     Null deviance: 208.15  on 158  degrees of freedom
## Residual deviance: 207.59  on 151  degrees of freedom
## AIC: 2794.7
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.4097 
##           Std. Err.:  0.0374 
## 
##  2 x log-likelihood:  -2776.6580 
## [1] "Model dispersion..."
## [1] 0.7875016
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
## 
## Response: n
##            LR Chisq Df Pr(>Chisq)
## modulation  0.55954  7     0.9992
## [1] "Incidence rates..."
##                   Estimate        2.5 %      97.5 %
## (Intercept)    2981.750000 1612.6030958 6456.973590
## modulation8Hz     1.043062    0.3885751    2.799918
## modulation16Hz    1.102859    0.4108522    2.960428
## modulation25Hz    1.130259    0.4210599    3.033976
## modulation33Hz    1.310958    0.4883784    3.519017
## modulation40Hz    1.214352    0.4482413    3.317702
## modulation50Hz    1.038920    0.3870321    2.788800
## modulation80Hz    1.260065    0.4694185    3.382408
## [1] "Estimated marginal mean contrasts..."
## $emmeans
##  modulation response   SE  df asymp.LCL asymp.UCL
##  4Hz            2982 1042 Inf      1504      5913
##  8Hz            3110 1087 Inf      1568      6168
##  16Hz           3288 1149 Inf      1658      6522
##  25Hz           3370 1177 Inf      1699      6684
##  33Hz           3909 1366 Inf      1971      7752
##  40Hz           3621 1298 Inf      1794      7310
##  50Hz           3098 1082 Inf      1562      6144
##  80Hz           3757 1313 Inf      1895      7451
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast    ratio    SE  df null z.ratio p.value
##  4Hz / 8Hz   0.959 0.474 Inf    1  -0.085  1.0000
##  4Hz / 16Hz  0.907 0.448 Inf    1  -0.198  1.0000
##  4Hz / 25Hz  0.885 0.437 Inf    1  -0.248  1.0000
##  4Hz / 33Hz  0.763 0.377 Inf    1  -0.548  1.0000
##  4Hz / 40Hz  0.823 0.412 Inf    1  -0.388  1.0000
##  4Hz / 50Hz  0.963 0.476 Inf    1  -0.077  1.0000
##  4Hz / 80Hz  0.794 0.392 Inf    1  -0.468  1.0000
##  8Hz / 16Hz  0.946 0.467 Inf    1  -0.113  1.0000
##  8Hz / 25Hz  0.923 0.456 Inf    1  -0.163  1.0000
##  8Hz / 33Hz  0.796 0.393 Inf    1  -0.463  1.0000
##  8Hz / 40Hz  0.859 0.430 Inf    1  -0.304  1.0000
##  8Hz / 50Hz  1.004 0.496 Inf    1   0.008  1.0000
##  8Hz / 80Hz  0.828 0.409 Inf    1  -0.383  1.0000
##  16Hz / 25Hz 0.976 0.482 Inf    1  -0.050  1.0000
##  16Hz / 33Hz 0.841 0.416 Inf    1  -0.350  1.0000
##  16Hz / 40Hz 0.908 0.455 Inf    1  -0.192  1.0000
##  16Hz / 50Hz 1.062 0.524 Inf    1   0.121  1.0000
##  16Hz / 80Hz 0.875 0.432 Inf    1  -0.270  1.0000
##  25Hz / 33Hz 0.862 0.426 Inf    1  -0.300  1.0000
##  25Hz / 40Hz 0.931 0.466 Inf    1  -0.143  1.0000
##  25Hz / 50Hz 1.088 0.537 Inf    1   0.171  1.0000
##  25Hz / 80Hz 0.897 0.443 Inf    1  -0.220  1.0000
##  33Hz / 40Hz 1.080 0.540 Inf    1   0.153  1.0000
##  33Hz / 50Hz 1.262 0.623 Inf    1   0.471  1.0000
##  33Hz / 80Hz 1.040 0.514 Inf    1   0.080  1.0000
##  40Hz / 50Hz 1.169 0.585 Inf    1   0.312  1.0000
##  40Hz / 80Hz 0.964 0.482 Inf    1  -0.074  1.0000
##  50Hz / 80Hz 0.824 0.407 Inf    1  -0.391  1.0000
## 
## P value adjustment: bonferroni method for 28 tests 
## Tests are performed on the log scale 
## 
## [1] "~~~~~~~~~~~~~~  all rates  ~~~~~~~~~~~~~~"
## [1] "Model summary..."
## 
## Call:
## MASS::glm.nb(formula = as.formula(formula), data = dat, init.theta = 0.3986165412, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3569  -1.4586  -0.6550   0.3964   1.7925  
## 
## Coefficients:
##                Estimate Std. Error z value            Pr(>|z|)    
## (Intercept)     7.27718    0.35422  20.545 <0.0000000000000002 ***
## modulation8Hz   0.86736    0.50746   1.709              0.0874 .  
## modulation16Hz  0.73719    0.50092   1.472              0.1411    
## modulation25Hz  0.17483    0.50748   0.345              0.7305    
## modulation33Hz  0.18570    0.50093   0.371              0.7109    
## modulation40Hz  0.21766    0.50748   0.429              0.6680    
## modulation50Hz -0.12942    0.50094  -0.258              0.7961    
## modulation80Hz  0.05829    0.50093   0.116              0.9074    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.3986) family taken to be 1)
## 
##     Null deviance: 212.84  on 156  degrees of freedom
## Residual deviance: 205.50  on 149  degrees of freedom
## AIC: 2566.6
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.3986 
##           Std. Err.:  0.0365 
## 
##  2 x log-likelihood:  -2548.6170 
## [1] "Model dispersion..."
## [1] 0.7789167
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
## 
## Response: n
##            LR Chisq Df Pr(>Chisq)
## modulation   7.3451  7     0.3939
## [1] "Incidence rates..."
##                    Estimate       2.5 %      97.5 %
## (Intercept)    1446.9000000 776.5128310 3171.476821
## modulation8Hz     2.3806250   0.8662196    6.599605
## modulation16Hz    2.0900546   0.7675921    5.690950
## modulation25Hz    1.1910400   0.4333624    3.301910
## modulation33Hz    1.2040569   0.4421913    3.278566
## modulation40Hz    1.2431660   0.4523296    3.446411
## modulation50Hz    0.8786025   0.3226615    2.392422
## modulation80Hz    1.0600249   0.3892925    2.886397
## [1] "Estimated marginal mean contrasts..."
## $emmeans
##  modulation response   SE  df asymp.LCL asymp.UCL
##  4Hz            1447  513 Inf       723      2897
##  8Hz            3445 1252 Inf      1690      7022
##  16Hz           3024 1071 Inf      1510      6055
##  25Hz           1723  626 Inf       845      3513
##  33Hz           1742  617 Inf       870      3488
##  40Hz           1799  654 Inf       882      3667
##  50Hz           1271  450 Inf       635      2545
##  80Hz           1534  543 Inf       766      3071
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast    ratio    SE  df null z.ratio p.value
##  4Hz / 8Hz   0.420 0.213 Inf    1  -1.709  1.0000
##  4Hz / 16Hz  0.478 0.240 Inf    1  -1.472  1.0000
##  4Hz / 25Hz  0.840 0.426 Inf    1  -0.345  1.0000
##  4Hz / 33Hz  0.831 0.416 Inf    1  -0.371  1.0000
##  4Hz / 40Hz  0.804 0.408 Inf    1  -0.429  1.0000
##  4Hz / 50Hz  1.138 0.570 Inf    1   0.258  1.0000
##  4Hz / 80Hz  0.943 0.473 Inf    1  -0.116  1.0000
##  8Hz / 16Hz  1.139 0.578 Inf    1   0.257  1.0000
##  8Hz / 25Hz  1.999 1.027 Inf    1   1.348  1.0000
##  8Hz / 33Hz  1.977 1.003 Inf    1   1.343  1.0000
##  8Hz / 40Hz  1.915 0.984 Inf    1   1.264  1.0000
##  8Hz / 50Hz  2.710 1.375 Inf    1   1.964  1.0000
##  8Hz / 80Hz  2.246 1.140 Inf    1   1.594  1.0000
##  16Hz / 25Hz 1.755 0.890 Inf    1   1.108  1.0000
##  16Hz / 33Hz 1.736 0.870 Inf    1   1.101  1.0000
##  16Hz / 40Hz 1.681 0.853 Inf    1   1.024  1.0000
##  16Hz / 50Hz 2.379 1.192 Inf    1   1.730  1.0000
##  16Hz / 80Hz 1.972 0.988 Inf    1   1.355  1.0000
##  25Hz / 33Hz 0.989 0.502 Inf    1  -0.021  1.0000
##  25Hz / 40Hz 0.958 0.492 Inf    1  -0.083  1.0000
##  25Hz / 50Hz 1.356 0.688 Inf    1   0.600  1.0000
##  25Hz / 80Hz 1.124 0.570 Inf    1   0.230  1.0000
##  33Hz / 40Hz 0.969 0.492 Inf    1  -0.063  1.0000
##  33Hz / 50Hz 1.370 0.686 Inf    1   0.629  1.0000
##  33Hz / 80Hz 1.136 0.569 Inf    1   0.254  1.0000
##  40Hz / 50Hz 1.415 0.718 Inf    1   0.684  1.0000
##  40Hz / 80Hz 1.173 0.595 Inf    1   0.314  1.0000
##  50Hz / 80Hz 0.829 0.415 Inf    1  -0.375  1.0000
## 
## P value adjustment: bonferroni method for 28 tests 
## Tests are performed on the log scale 
## 
## [1] "~~~~~~~~~~~~~~  all rates  ~~~~~~~~~~~~~~"
## [1] "Model summary..."
## 
## Call:
## MASS::glm.nb(formula = as.formula(formula), data = dat, init.theta = 0.5412110414, 
##     link = log)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.91889  -1.35093  -0.49683   0.08802   1.60111  
## 
## Coefficients:
##                Estimate Std. Error z value             Pr(>|z|)    
## (Intercept)      8.5754     0.3040  28.212 < 0.0000000000000002 ***
## modulation8Hz   -0.7133     0.4299  -1.659             0.097070 .  
## modulation16Hz  -1.4230     0.4299  -3.310             0.000933 ***
## modulation25Hz  -1.8539     0.4299  -4.312     0.00001616474449 ***
## modulation33Hz  -2.1116     0.4300  -4.911     0.00000090540685 ***
## modulation40Hz  -2.3090     0.4300  -5.370     0.00000007865917 ***
## modulation50Hz  -2.5248     0.4300  -5.872     0.00000000431401 ***
## modulation80Hz  -3.0032     0.4301  -6.983     0.00000000000289 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.5412) family taken to be 1)
## 
##     Null deviance: 281.09  on 159  degrees of freedom
## Residual deviance: 200.02  on 152  degrees of freedom
## AIC: 2473.2
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.5412 
##           Std. Err.:  0.0505 
## 
##  2 x log-likelihood:  -2455.2150 
## [1] "Model dispersion..."
## [1] 1.06478
## [1] "Type II Deviance Table..."
## Analysis of Deviance Table (Type II tests)
## 
## Response: n
##            LR Chisq Df         Pr(>Chisq)    
## modulation   81.078  7 0.0000000000000083 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Incidence rates..."
##                     Estimate         2.5 %        97.5 %
## (Intercept)    5299.50000000 3082.24798764 10267.1812755
## modulation8Hz     0.49003680    0.20837906     1.1524002
## modulation16Hz    0.24099443    0.10247397     0.5667617
## modulation25Hz    0.15661855    0.06659318     0.3683466
## modulation33Hz    0.12104916    0.05146732     0.2847030
## modulation40Hz    0.09935843    0.04224335     0.2336959
## modulation50Hz    0.08007359    0.03404249     0.1883464
## modulation80Hz    0.04962732    0.02109523     0.1167501
## [1] "Estimated marginal mean contrasts..."
## $emmeans
##  modulation response   SE  df asymp.LCL asymp.UCL
##  4Hz            5300 1611 Inf      2921      9615
##  8Hz            2597  789 Inf      1431      4712
##  16Hz           1277  388 Inf       704      2317
##  25Hz            830  252 Inf       457      1506
##  33Hz            642  195 Inf       353      1164
##  40Hz            527  160 Inf       290       956
##  50Hz            424  129 Inf       234       770
##  80Hz            263   80 Inf       145       477
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale 
## 
## $contrasts
##  contrast    ratio    SE  df null z.ratio p.value
##  4Hz / 8Hz    2.04 0.877 Inf    1   1.659  1.0000
##  4Hz / 16Hz   4.15 1.784 Inf    1   3.310  0.0261
##  4Hz / 25Hz   6.38 2.745 Inf    1   4.312  0.0005
##  4Hz / 33Hz   8.26 3.552 Inf    1   4.911  <.0001
##  4Hz / 40Hz  10.06 4.327 Inf    1   5.370  <.0001
##  4Hz / 50Hz  12.49 5.370 Inf    1   5.872  <.0001
##  4Hz / 80Hz  20.15 8.666 Inf    1   6.983  <.0001
##  8Hz / 16Hz   2.03 0.874 Inf    1   1.651  1.0000
##  8Hz / 25Hz   3.13 1.345 Inf    1   2.653  0.2233
##  8Hz / 33Hz   4.05 1.741 Inf    1   3.252  0.0321
##  8Hz / 40Hz   4.93 2.121 Inf    1   3.711  0.0058
##  8Hz / 50Hz   6.12 2.632 Inf    1   4.213  0.0007
##  8Hz / 80Hz   9.87 4.247 Inf    1   5.324  <.0001
##  16Hz / 25Hz  1.54 0.662 Inf    1   1.002  1.0000
##  16Hz / 33Hz  1.99 0.856 Inf    1   1.601  1.0000
##  16Hz / 40Hz  2.43 1.043 Inf    1   2.061  1.0000
##  16Hz / 50Hz  3.01 1.294 Inf    1   2.562  0.2912
##  16Hz / 80Hz  4.86 2.089 Inf    1   3.674  0.0067
##  25Hz / 33Hz  1.29 0.556 Inf    1   0.599  1.0000
##  25Hz / 40Hz  1.58 0.678 Inf    1   1.058  1.0000
##  25Hz / 50Hz  1.96 0.841 Inf    1   1.560  1.0000
##  25Hz / 80Hz  3.16 1.357 Inf    1   2.672  0.2112
##  33Hz / 40Hz  1.22 0.524 Inf    1   0.459  1.0000
##  33Hz / 50Hz  1.51 0.650 Inf    1   0.961  1.0000
##  33Hz / 80Hz  2.44 1.049 Inf    1   2.073  1.0000
##  40Hz / 50Hz  1.24 0.534 Inf    1   0.502  1.0000
##  40Hz / 80Hz  2.00 0.861 Inf    1   1.614  1.0000
##  50Hz / 80Hz  1.61 0.694 Inf    1   1.112  1.0000
## 
## P value adjustment: bonferroni method for 28 tests 
## Tests are performed on the log scale

names(all.models2) <- sort(unique(n_calls_reg$condition))

Anova

# pull out anova tables
(anova2 <- do.call(rbind, lapply(all.models2, '[[', 'deviance')) %>% 
   mutate(across(is.numeric, round, 2)))

##                     LR Chisq Df Pr(>Chisq)
## silence                 0.56  7       1.00
## steady-state masker     7.34  7       0.39
## random masker          81.08  7       0.00

Marginal mean contrasts

##             modulation    contrast   IRR   se asymp_lcl asymp_ucl z_ratio p_value
## 1              silence   4Hz / 8Hz  0.96 0.47      0.20      4.49   -0.09    1.00
## 2              silence  4Hz / 16Hz  0.91 0.45      0.19      4.24   -0.20    1.00
## 3              silence  4Hz / 25Hz  0.88 0.44      0.19      4.14   -0.25    1.00
## 4              silence  4Hz / 33Hz  0.76 0.38      0.16      3.57   -0.55    1.00
## 5              silence  4Hz / 40Hz  0.82 0.41      0.17      3.93   -0.39    1.00
## 6              silence  4Hz / 50Hz  0.96 0.48      0.21      4.50   -0.08    1.00
## 7              silence  4Hz / 80Hz  0.79 0.39      0.17      3.71   -0.47    1.00
## 8              silence  8Hz / 16Hz  0.95 0.47      0.20      4.43   -0.11    1.00
## 9              silence  8Hz / 25Hz  0.92 0.46      0.20      4.32   -0.16    1.00
## 10             silence  8Hz / 33Hz  0.80 0.39      0.17      3.72   -0.46    1.00
## 11             silence  8Hz / 40Hz  0.86 0.43      0.18      4.10   -0.30    1.00
## 12             silence  8Hz / 50Hz  1.00 0.50      0.22      4.70    0.01    1.00
## 13             silence  8Hz / 80Hz  0.83 0.41      0.18      3.87   -0.38    1.00
## 14             silence 16Hz / 25Hz  0.98 0.48      0.21      4.57   -0.05    1.00
## 15             silence 16Hz / 33Hz  0.84 0.42      0.18      3.94   -0.35    1.00
## 16             silence 16Hz / 40Hz  0.91 0.46      0.19      4.34   -0.19    1.00
## 17             silence 16Hz / 50Hz  1.06 0.52      0.23      4.97    0.12    1.00
## 18             silence 16Hz / 80Hz  0.88 0.43      0.19      4.10   -0.27    1.00
## 19             silence 25Hz / 33Hz  0.86 0.43      0.18      4.04   -0.30    1.00
## 20             silence 25Hz / 40Hz  0.93 0.47      0.20      4.44   -0.14    1.00
## 21             silence 25Hz / 50Hz  1.09 0.54      0.23      5.09    0.17    1.00
## 22             silence 25Hz / 80Hz  0.90 0.44      0.19      4.20   -0.22    1.00
## 23             silence 33Hz / 40Hz  1.08 0.54      0.23      5.16    0.15    1.00
## 24             silence 33Hz / 50Hz  1.26 0.62      0.27      5.90    0.47    1.00
## 25             silence 33Hz / 80Hz  1.04 0.51      0.22      4.87    0.08    1.00
## 26             silence 40Hz / 50Hz  1.17 0.58      0.24      5.58    0.31    1.00
## 27             silence 40Hz / 80Hz  0.96 0.48      0.20      4.60   -0.07    1.00
## 28             silence 50Hz / 80Hz  0.82 0.41      0.18      3.86   -0.39    1.00
## 29 steady-state masker   4Hz / 8Hz  0.42 0.21      0.09      2.05   -1.71    1.00
## 30 steady-state masker  4Hz / 16Hz  0.48 0.24      0.10      2.29   -1.47    1.00
## 31 steady-state masker  4Hz / 25Hz  0.84 0.43      0.17      4.10   -0.34    1.00
## 32 steady-state masker  4Hz / 33Hz  0.83 0.42      0.17      3.97   -0.37    1.00
## 33 steady-state masker  4Hz / 40Hz  0.80 0.41      0.16      3.93   -0.43    1.00
## 34 steady-state masker  4Hz / 50Hz  1.14 0.57      0.24      5.44    0.26    1.00
## 35 steady-state masker  4Hz / 80Hz  0.94 0.47      0.20      4.51   -0.12    1.00
## 36 steady-state masker  8Hz / 16Hz  1.14 0.58      0.23      5.56    0.26    1.00
## 37 steady-state masker  8Hz / 25Hz  2.00 1.03      0.40      9.95    1.35    1.00
## 38 steady-state masker  8Hz / 33Hz  1.98 1.00      0.41      9.65    1.34    1.00
## 39 steady-state masker  8Hz / 40Hz  1.92 0.98      0.38      9.54    1.26    1.00
## 40 steady-state masker  8Hz / 50Hz  2.71 1.38      0.56     13.22    1.96    1.00
## 41 steady-state masker  8Hz / 80Hz  2.25 1.14      0.46     10.96    1.59    1.00
## 42 steady-state masker 16Hz / 25Hz  1.75 0.89      0.36      8.56    1.11    1.00
## 43 steady-state masker 16Hz / 33Hz  1.74 0.87      0.36      8.30    1.10    1.00
## 44 steady-state masker 16Hz / 40Hz  1.68 0.85      0.34      8.20    1.02    1.00
## 45 steady-state masker 16Hz / 50Hz  2.38 1.19      0.50     11.37    1.73    1.00
## 46 steady-state masker 16Hz / 80Hz  1.97 0.99      0.41      9.43    1.35    1.00
## 47 steady-state masker 25Hz / 33Hz  0.99 0.50      0.20      4.83   -0.02    1.00
## 48 steady-state masker 25Hz / 40Hz  0.96 0.49      0.19      4.77   -0.08    1.00
## 49 steady-state masker 25Hz / 50Hz  1.36 0.69      0.28      6.62    0.60    1.00
## 50 steady-state masker 25Hz / 80Hz  1.12 0.57      0.23      5.48    0.23    1.00
## 51 steady-state masker 33Hz / 40Hz  0.97 0.49      0.20      4.73   -0.06    1.00
## 52 steady-state masker 33Hz / 50Hz  1.37 0.69      0.29      6.55    0.63    1.00
## 53 steady-state masker 33Hz / 80Hz  1.14 0.57      0.24      5.43    0.25    1.00
## 54 steady-state masker 40Hz / 50Hz  1.42 0.72      0.29      6.91    0.68    1.00
## 55 steady-state masker 40Hz / 80Hz  1.17 0.60      0.24      5.72    0.31    1.00
## 56 steady-state masker 50Hz / 80Hz  0.83 0.42      0.17      3.96   -0.38    1.00
## 57       random masker   4Hz / 8Hz  2.04 0.88      0.53      7.82    1.66    1.00
## 58       random masker  4Hz / 16Hz  4.15 1.78      1.08     15.89    3.31    0.03
## 59       random masker  4Hz / 25Hz  6.38 2.74      1.67     24.46    4.31    0.00
## 60       random masker  4Hz / 33Hz  8.26 3.55      2.16     31.65    4.91    0.00
## 61       random masker  4Hz / 40Hz 10.06 4.33      2.63     38.56    5.37    0.00
## 62       random masker  4Hz / 50Hz 12.49 5.37      3.26     47.85    5.87    0.00
## 63       random masker  4Hz / 80Hz 20.15 8.67      5.26     77.22    6.98    0.00
## 64       random masker  8Hz / 16Hz  2.03 0.87      0.53      7.79    1.65    1.00
## 65       random masker  8Hz / 25Hz  3.13 1.34      0.82     11.98    2.65    0.22
## 66       random masker  8Hz / 33Hz  4.05 1.74      1.06     15.51    3.25    0.03
## 67       random masker  8Hz / 40Hz  4.93 2.12      1.29     18.90    3.71    0.01
## 68       random masker  8Hz / 50Hz  6.12 2.63      1.60     23.45    4.21    0.00
## 69       random masker  8Hz / 80Hz  9.87 4.25      2.58     37.84    5.32    0.00
## 70       random masker 16Hz / 25Hz  1.54 0.66      0.40      5.89    1.00    1.00
## 71       random masker 16Hz / 33Hz  1.99 0.86      0.52      7.63    1.60    1.00
## 72       random masker 16Hz / 40Hz  2.43 1.04      0.63      9.29    2.06    1.00
## 73       random masker 16Hz / 50Hz  3.01 1.29      0.78     11.53    2.56    0.29
## 74       random masker 16Hz / 80Hz  4.86 2.09      1.27     18.61    3.67    0.01
## 75       random masker 25Hz / 33Hz  1.29 0.56      0.34      4.96    0.60    1.00
## 76       random masker 25Hz / 40Hz  1.58 0.68      0.41      6.04    1.06    1.00
## 77       random masker 25Hz / 50Hz  1.96 0.84      0.51      7.50    1.56    1.00
## 78       random masker 25Hz / 80Hz  3.16 1.36      0.82     12.10    2.67    0.21
## 79       random masker 33Hz / 40Hz  1.22 0.52      0.32      4.67    0.46    1.00
## 80       random masker 33Hz / 50Hz  1.51 0.65      0.39      5.79    0.96    1.00
## 81       random masker 33Hz / 80Hz  2.44 1.05      0.64      9.35    2.07    1.00
## 82       random masker 40Hz / 50Hz  1.24 0.53      0.32      4.76    0.50    1.00
## 83       random masker 40Hz / 80Hz  2.00 0.86      0.52      7.68    1.61    1.00
## 84       random masker 50Hz / 80Hz  1.61 0.69      0.42      6.19    1.11    1.00

Incidence Rate Ratios (IRR)

# IRR
IRRs2 <- do.call(rbind, lapply(all.models2, '[[', 'IRR')) %>%
   mutate(across(is.numeric, round, 2)) %>%
   mutate(condition = str_split(row.names(.),"[.]",simplify = T)[,1], .before=1) %>%
   mutate(coef = str_remove(str_split(row.names(.),"[.]",simplify = T)[,2],"condition"), .after=1) %>%
   remove_rownames()
IRRs2

##              condition           coef Estimate   2.5 %   97.5 %
## 1              silence    (Intercept)  2981.75 1612.60  6456.97
## 2              silence  modulation8Hz     1.04    0.39     2.80
## 3              silence modulation16Hz     1.10    0.41     2.96
## 4              silence modulation25Hz     1.13    0.42     3.03
## 5              silence modulation33Hz     1.31    0.49     3.52
## 6              silence modulation40Hz     1.21    0.45     3.32
## 7              silence modulation50Hz     1.04    0.39     2.79
## 8              silence modulation80Hz     1.26    0.47     3.38
## 9  steady-state masker    (Intercept)  1446.90  776.51  3171.48
## 10 steady-state masker  modulation8Hz     2.38    0.87     6.60
## 11 steady-state masker modulation16Hz     2.09    0.77     5.69
## 12 steady-state masker modulation25Hz     1.19    0.43     3.30
## 13 steady-state masker modulation33Hz     1.20    0.44     3.28
## 14 steady-state masker modulation40Hz     1.24    0.45     3.45
## 15 steady-state masker modulation50Hz     0.88    0.32     2.39
## 16 steady-state masker modulation80Hz     1.06    0.39     2.89
## 17       random masker    (Intercept)  5299.50 3082.25 10267.18
## 18       random masker  modulation8Hz     0.49    0.21     1.15
## 19       random masker modulation16Hz     0.24    0.10     0.57
## 20       random masker modulation25Hz     0.16    0.07     0.37
## 21       random masker modulation33Hz     0.12    0.05     0.28
## 22       random masker modulation40Hz     0.10    0.04     0.23
## 23       random masker modulation50Hz     0.08    0.03     0.19
## 24       random masker modulation80Hz     0.05    0.02     0.12

Anticipating the trough or the peak?

Use bootstrapped means from the whole data set, we now want to see whether average call onsets were “aimed” more at dodging AM peaks or more targeted at AM troughs.

To find out, we first compute two measures of call onset timing as follows:

Plot timings relative to both acoustic features. Show raw data (dots), distributions, and boxplots.

# timings by condition and modulation
more_timings <- mu_bs_time %>% group_by(modulation, condition) %>% 
  summarise(median_peak = median(t_from_peak), 
           iqr_peak = IQR(t_from_peak),
           mad_peak = mad(t_from_peak),
           median_trough = median(t_to_trough),
           mad_trough = mad(t_to_trough),
           iqr_trough = IQR(t_to_trough))

more_timings

## # A tibble: 24 × 8
## # Groups:   modulation [8]
##    modulation condition           median_peak iqr_peak mad_peak median_trough mad_trough iqr_trough
##    <fct>      <fct>                     <dbl>    <dbl>    <dbl>         <dbl>      <dbl>      <dbl>
##  1 4Hz        silence                  0.0308 0.200    0.131           0.0942   0.131      0.200   
##  2 4Hz        steady-state masker      0.0837 0.00310  0.00222         0.0413   0.00222    0.00310 
##  3 4Hz        random masker            0.0837 0.00210  0.00163         0.0413   0.00163    0.00210 
##  4 8Hz        silence                  0.0346 0.0456   0.0249          0.0279   0.0249     0.0456  
##  5 8Hz        steady-state masker      0.0528 0.00130  0.000890        0.0097   0.000890   0.00130 
##  6 8Hz        random masker            0.0462 0.00140  0.00104         0.0163   0.00104    0.00140 
##  7 16Hz       silence                 -0.0197 0.0144   0.00875         0.0505   0.00875    0.0144  
##  8 16Hz       steady-state masker      0.0258 0.000900 0.000593        0.0050   0.000593   0.000900
##  9 16Hz       random masker            0.0223 0.00140  0.00104         0.0085   0.00104    0.00140 
## 10 25Hz       silence                  0.0099 0.00663  0.00482         0.0101   0.00482    0.00663 
## # … with 14 more rows

pb <- p.peak_trough(mu_bs_time %>% 
                      dplyr::filter(!condition=="silence", !condition=="half-band masker"), binwidth = 8) %>%
  recolor()

print(pb)

ggsave(file.path(figures_path,paste0(exp, '_fig3a-alt1.png')), pb, width = 6, height = 4)

Random Forest models

Run three versions of the random forest model:

1. A “full” model: \(modulation ~ t_{peak} + t_{trough}\)

2. A “peaks” model: \(modulation ~ t_{peak}\)

3. A “troughs” model: \(modulation ~ t_{trough}\)

mu_rf <- mu_bs_time %>% dplyr::filter(!condition=="silence") %>%
  dplyr::select(t_to_trough, t_from_peak, modulation)

test_idx <- sample(1:nrow(mu_rf), size = floor(nrow(mu_rf)*0.60))

# split into 0.6:0.4 train/validation set
mu_bs_train <- mu_rf[test_idx,] 
mu_bs_valid <- mu_rf[-test_idx,] 

# shuffled
mu_bs_shuffled_pk <- mu_bs_train

rf_files <- c(file.path(models_path, paste0(exp,"_RF_full.RDS")),
              file.path(models_path, paste0(exp,"_RF_troughs_only.RDS")),
              file.path(models_path, paste0(exp,"_RF_peaks_only.RDS")))



if (all(file.exists(rf_files))) {
  rf_models <- lapply(rf_files, FUN = readr::read_rds)
  rf_classifier <- rf_models[[1]]
  rf_classifier_tr <- rf_models[[2]]
  rf_classifier_pk <- rf_models[[3]]
} else {

  # train
  # full model
  rf_classifier <- randomForest(modulation ~ ., data=mu_bs_train, ntree=500, mtry=2, importance=TRUE, proximity=TRUE)
  readr::write_rds(rf_classifier, file = rf_files[1])
  
  # only troughs
  rf_classifier_tr <- randomForest(modulation ~ t_to_trough, data=mu_bs_train, ntree=500, mtry=1, proximity=TRUE)
  readr::write_rds(rf_classifier_tr, file = rf_files[2])
  
  # only peaks
  rf_classifier_pk <- randomForest(modulation ~ t_from_peak, data=mu_bs_train, ntree=500, mtry=1, proximity=TRUE)
  readr::write_rds(rf_classifier_pk, file = rf_files[3])
  
  rf_models <- list(rf_classifier, rf_classifier_tr, rf_classifier_pk)


}
rf_classifier

## 
## Call:
##  randomForest(formula = modulation ~ ., data = mu_bs_train, ntree = 500,      mtry = 2, importance = TRUE, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##         OOB estimate of  error rate: 0.26%
## Confusion matrix:
##       4Hz  8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz class.error
## 4Hz  1183    0    0    0    0    0    0    0 0.000000000
## 8Hz     0 1245    0    0    0    0    0    0 0.000000000
## 16Hz    0    0 1179    0    0    0    0    0 0.000000000
## 25Hz    0    0    0 1192   15    0    0    0 0.012427506
## 33Hz    0    0    0    0 1184    2    0    0 0.001686341
## 40Hz    0    0    0    0    2 1248    3    0 0.003990423
## 50Hz    0    0    0    0    0    0 1154    3 0.002592913
## 80Hz    0    0    0    0    0    0    0 1190 0.000000000

rf_classifier_tr

## 
## Call:
##  randomForest(formula = modulation ~ t_to_trough, data = mu_bs_train,      ntree = 500, mtry = 1, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 1
## 
##         OOB estimate of  error rate: 36.68%
## Confusion matrix:
##       4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz class.error
## 4Hz  1183   0    0    0    0    0    0    0   0.0000000
## 8Hz     0 848   40   16    5   80  229   27   0.3188755
## 16Hz    0  52  760    8   73  128  108   50   0.3553859
## 25Hz    0   9   99  280   86  174  342  217   0.7680199
## 33Hz    0   6  212    4  697   59   20  188   0.4123103
## 40Hz    0   1   58    1   37  652  223  281   0.4796488
## 50Hz    0  70   17   63   17   36  744  210   0.3569576
## 80Hz    0  11   86   38   79   14   47  915   0.2310924

rf_classifier_pk

## 
## Call:
##  randomForest(formula = modulation ~ t_from_peak, data = mu_bs_train,      ntree = 500, mtry = 1, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 1
## 
##         OOB estimate of  error rate: 15.11%
## Confusion matrix:
##       4Hz  8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz class.error
## 4Hz  1183    0    0    0    0    0    0    0  0.00000000
## 8Hz     0 1245    0    0    0    0    0    0  0.00000000
## 16Hz    0    0 1179    0    0    0    0    0  0.00000000
## 25Hz    0    0    0  730  436    5   36    0  0.39519470
## 33Hz    0    0    0    5 1129    2    2   48  0.04806071
## 40Hz    0    0    0   34   13 1010   49  147  0.19393456
## 50Hz    0    0    0    1   32   18  994  112  0.14088159
## 80Hz    0    0    0    0   12  202  297  679  0.42941176

Calculate variable importance of two predictors from the full model…

Predict modulation classes for the validation set for each of the models…

predictions <- predict(rf_classifier,mu_bs_valid[,-3])
predictionstr <- predict(rf_classifier_tr,mu_bs_valid[,-3])
predictionspk <- predict(rf_classifier_pk,mu_bs_valid[,-3])

# predicted classes vs. observed classes
table(observed=mu_bs_valid[,3],predicted=predictions)

##         predicted
## observed 4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz
##     4Hz  817   0    0    0    0    0    0    0
##     8Hz    0 755    0    0    0    0    0    0
##     16Hz   0   0  821    0    0    0    0    0
##     25Hz   0   0    0  783   10    0    0    0
##     33Hz   0   0    0    0  811    3    0    0
##     40Hz   0   0    0    0    2  745    0    0
##     50Hz   0   0    0    0    0    2  838    3
##     80Hz   0   0    0    0    0    0    0  810

print("Full model:")

## [1] "Full model:"

# proportion misclassifications
mean(predictions != mu_bs_valid$modulation) * 100

## [1] 0.3125

# stats
(matfull <- caret::confusionMatrix(predictions, mu_bs_valid$modulation))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction 4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz
##       4Hz  817   0    0    0    0    0    0    0
##       8Hz    0 755    0    0    0    0    0    0
##       16Hz   0   0  821    0    0    0    0    0
##       25Hz   0   0    0  783    0    0    0    0
##       33Hz   0   0    0   10  811    2    0    0
##       40Hz   0   0    0    0    3  745    2    0
##       50Hz   0   0    0    0    0    0  838    0
##       80Hz   0   0    0    0    0    0    3  810
## 
## Overall Statistics
##                                                
##                Accuracy : 0.9969               
##                  95% CI : (0.9952, 0.9981)     
##     No Information Rate : 0.1317               
##     P-Value [Acc > NIR] : < 0.00000000000000022
##                                                
##                   Kappa : 0.9964               
##                                                
##  Mcnemar's Test P-Value : NA                   
## 
## Statistics by Class:
## 
##                      Class: 4Hz Class: 8Hz Class: 16Hz Class: 25Hz Class: 33Hz Class: 40Hz Class: 50Hz
## Sensitivity              1.0000      1.000      1.0000      0.9874      0.9963      0.9973      0.9941
## Specificity              1.0000      1.000      1.0000      1.0000      0.9979      0.9991      1.0000
## Pos Pred Value           1.0000      1.000      1.0000      1.0000      0.9854      0.9933      1.0000
## Neg Pred Value           1.0000      1.000      1.0000      0.9982      0.9995      0.9996      0.9991
## Prevalence               0.1277      0.118      0.1283      0.1239      0.1272      0.1167      0.1317
## Detection Rate           0.1277      0.118      0.1283      0.1223      0.1267      0.1164      0.1309
## Detection Prevalence     0.1277      0.118      0.1283      0.1223      0.1286      0.1172      0.1309
## Balanced Accuracy        1.0000      1.000      1.0000      0.9937      0.9971      0.9982      0.9970
##                      Class: 80Hz
## Sensitivity               1.0000
## Specificity               0.9995
## Pos Pred Value            0.9963
## Neg Pred Value            1.0000
## Prevalence                0.1266
## Detection Rate            0.1266
## Detection Prevalence      0.1270
## Balanced Accuracy         0.9997

print("Troughs model:")

## [1] "Troughs model:"

# stats
(mattr <- caret::confusionMatrix(predictionstr, mu_bs_valid$modulation))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction 4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz
##       4Hz  817   0    0    0    0    0    0    0
##       8Hz    0 520   35    6    4    0   37   12
##       16Hz   0  23  539   64  156   39   14   41
##       25Hz   0  13    3  193    4    1   52   39
##       33Hz   0   3   56   50  449    9   15   65
##       40Hz   0  49   91  104   49  390   26    5
##       50Hz   0 123   67  215   17  130  568   41
##       80Hz   0  24   30  161  135  178  131  607
## 
## Overall Statistics
##                                                
##                Accuracy : 0.638                
##                  95% CI : (0.6261, 0.6498)     
##     No Information Rate : 0.1317               
##     P-Value [Acc > NIR] : < 0.00000000000000022
##                                                
##                   Kappa : 0.5858               
##                                                
##  Mcnemar's Test P-Value : NA                   
## 
## Statistics by Class:
## 
##                      Class: 4Hz Class: 8Hz Class: 16Hz Class: 25Hz Class: 33Hz Class: 40Hz Class: 50Hz
## Sensitivity              1.0000    0.68874     0.65652     0.24338     0.55160     0.52209     0.67378
## Specificity              1.0000    0.98335     0.93959     0.98002     0.96455     0.94269     0.89329
## Pos Pred Value           1.0000    0.84691     0.61530     0.63279     0.69397     0.54622     0.48923
## Neg Pred Value           1.0000    0.95938     0.94895     0.90156     0.93655     0.93721     0.94751
## Prevalence               0.1277    0.11797     0.12828     0.12391     0.12719     0.11672     0.13172
## Detection Rate           0.1277    0.08125     0.08422     0.03016     0.07016     0.06094     0.08875
## Detection Prevalence     0.1277    0.09594     0.13687     0.04766     0.10109     0.11156     0.18141
## Balanced Accuracy        1.0000    0.83604     0.79806     0.61170     0.75808     0.73239     0.78354
##                      Class: 80Hz
## Sensitivity              0.74938
## Specificity              0.88211
## Pos Pred Value           0.47946
## Neg Pred Value           0.96046
## Prevalence               0.12656
## Detection Rate           0.09484
## Detection Prevalence     0.19781
## Balanced Accuracy        0.81575

print("Peaks model:")

## [1] "Peaks model:"

# stats
(matpk <- caret::confusionMatrix(predictionspk, mu_bs_valid$modulation))

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction 4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz
##       4Hz  817   0    0    0    0    0    0    0
##       8Hz    0 755    0    0    0    0    0    0
##       16Hz   0   0  821    0    0    0    0    0
##       25Hz   0   0    0  471    4   26    0    0
##       33Hz   0   0    0  290  771    8   19    0
##       40Hz   0   0    0    8    0  603   17  150
##       50Hz   0   0    0   23    4   33  747  190
##       80Hz   0   0    0    1   35   77   60  470
## 
## Overall Statistics
##                                                
##                Accuracy : 0.8523               
##                  95% CI : (0.8434, 0.861)      
##     No Information Rate : 0.1317               
##     P-Value [Acc > NIR] : < 0.00000000000000022
##                                                
##                   Kappa : 0.8312               
##                                                
##  Mcnemar's Test P-Value : NA                   
## 
## Statistics by Class:
## 
##                      Class: 4Hz Class: 8Hz Class: 16Hz Class: 25Hz Class: 33Hz Class: 40Hz Class: 50Hz
## Sensitivity              1.0000      1.000      1.0000     0.59395      0.9472     0.80723      0.8861
## Specificity              1.0000      1.000      1.0000     0.99465      0.9433     0.96904      0.9550
## Pos Pred Value           1.0000      1.000      1.0000     0.94012      0.7086     0.77506      0.7492
## Neg Pred Value           1.0000      1.000      1.0000     0.94541      0.9919     0.97439      0.9822
## Prevalence               0.1277      0.118      0.1283     0.12391      0.1272     0.11672      0.1317
## Detection Rate           0.1277      0.118      0.1283     0.07359      0.1205     0.09422      0.1167
## Detection Prevalence     0.1277      0.118      0.1283     0.07828      0.1700     0.12156      0.1558
## Balanced Accuracy        1.0000      1.000      1.0000     0.79430      0.9452     0.88814      0.9206
##                      Class: 80Hz
## Sensitivity              0.58025
## Specificity              0.96905
## Pos Pred Value           0.73095
## Neg Pred Value           0.94094
## Prevalence               0.12656
## Detection Rate           0.07344
## Detection Prevalence     0.10047
## Balanced Accuracy        0.77465

validation <- list(matfull, mattr, matpk)

Plot confusion matrices from training models

# training matrices
# confmats <- lapply(rf_models,function(x) cormat_to_df(x$confusion[,-3])) %>%
#   setNames(c("full","troughs","peaks")) %>%
#   map_df(., ~as.data.frame(.x), .id="model") %>%
#   mutate(model = factor(model, levels = c("troughs","peaks","full")))

# validation matrices
confmats <- lapply(validation, function(x) cormat_to_df(x$table %>% as.data.frame()) ) %>%
  setNames(c("full", "troughs", "peaks")) %>%
  map_df(., ~as.data.frame(.x), .id="model") %>%
  mutate(model = factor(model, levels = c("troughs", "peaks", "full")))


confmat <- p.mat(confmats, ulim = max((mu_bs_valid%>%group_by(modulation)%>%summarise(n=n()))[,2])) +
  facet_grid(cols=vars(model)) 
  print(confmat)

  ggsave(file.path(figures_path,paste0(exp, '_fig3b.png')), 
         confmat, width = h, height = w2) 

Shuffled Forests

Now, let’s see what happens if we shuffle the modulation class labelled for one or the other measure in our full model. The drop in performance, respectively, should give us some insight into which variable is doing more of the legwork.

# shuffled
mu_bs_shuffled_pk <- mu_bs_train
mu_bs_shuffled_pk$t_from_peak <- mu_bs_shuffled_pk$t_from_peak[sample(1:nrow(mu_bs_shuffled_pk))]

mu_bs_shuffled_tr <- mu_bs_train
mu_bs_shuffled_tr$t_to_trough <- mu_bs_shuffled_tr$t_to_trough[sample(1:nrow(mu_bs_shuffled_tr))]

# train
# full model
shPeak_rf_classifier <- randomForest(modulation ~ ., data=mu_bs_shuffled_pk, ntree=500, mtry=2, importance=TRUE, proximity=TRUE)
shPeak_rf_classifier

## 
## Call:
##  randomForest(formula = modulation ~ ., data = mu_bs_shuffled_pk,      ntree = 500, mtry = 2, importance = TRUE, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##         OOB estimate of  error rate: 42.24%
## Confusion matrix:
##       4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz class.error
## 4Hz  1183   0    0    0    0    0    0    0   0.0000000
## 8Hz     0 865   58   21   13   57  197   34   0.3052209
## 16Hz    0  60  569   91  171   82  108   98   0.5173876
## 25Hz    0  19   98  299   76  112  328  275   0.7522784
## 33Hz    0   9  162   81  601   33   32  268   0.4932546
## 40Hz    0   5   93  116   37  489  210  303   0.6097366
## 50Hz    0 105   17   69   30   42  655  239   0.4338807
## 80Hz    0   8   49   48  131   13   57  884   0.2571429

shTrough_rf_classifier <- randomForest(modulation ~ ., data=mu_bs_shuffled_tr, ntree=500, mtry=2, importance=TRUE, proximity=TRUE)
shTrough_rf_classifier

## 
## Call:
##  randomForest(formula = modulation ~ ., data = mu_bs_shuffled_tr,      ntree = 500, mtry = 2, importance = TRUE, proximity = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##         OOB estimate of  error rate: 17.32%
## Confusion matrix:
##       4Hz  8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz class.error
## 4Hz  1183    0    0    0    0    0    0    0   0.0000000
## 8Hz     0 1245    0    0    0    0    0    0   0.0000000
## 16Hz    0    0 1179    0    0    0    0    0   0.0000000
## 25Hz    0    0    0  757  411   30    9    0   0.3728252
## 33Hz    0    0    0   18 1095    6   29   38   0.0767285
## 40Hz    0    0    0   26   14  962   52  199   0.2322426
## 50Hz    0    0    0    1   34   35  882  205   0.2376837
## 80Hz    0    0    0    0   33  224  299  634   0.4672269

Let’s get predictions as before…

predictionstr_sh <- predict(shTrough_rf_classifier,mu_bs_valid[,-3])
predictionspk_sh <- predict(shPeak_rf_classifier,mu_bs_valid[,-3])

print("Shuffled troughs:")

## [1] "Shuffled troughs:"

# stats
shtr <- caret::confusionMatrix(predictionstr_sh , mu_bs_valid$modulation)
shtr

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction 4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz
##       4Hz  817   0    0    0    0    0    0    0
##       8Hz    0 755    0    0    0    0    0    0
##       16Hz   0   0  821    0    0    0    0    0
##       25Hz   0   0    0  502   24   24    1    0
##       33Hz   0   0    0  264  744    7   19   28
##       40Hz   0   0    0    8    2  616   27  111
##       50Hz   0   0    0   18   22   23  637  178
##       80Hz   0   0    0    1   22   77  159  493
## 
## Overall Statistics
##                                                
##                Accuracy : 0.8414               
##                  95% CI : (0.8322, 0.8503)     
##     No Information Rate : 0.1317               
##     P-Value [Acc > NIR] : < 0.00000000000000022
##                                                
##                   Kappa : 0.8187               
##                                                
##  Mcnemar's Test P-Value : NA                   
## 
## Statistics by Class:
## 
##                      Class: 4Hz Class: 8Hz Class: 16Hz Class: 25Hz Class: 33Hz Class: 40Hz Class: 50Hz
## Sensitivity              1.0000      1.000      1.0000     0.63304      0.9140     0.82463     0.75563
## Specificity              1.0000      1.000      1.0000     0.99126      0.9431     0.97382     0.95663
## Pos Pred Value           1.0000      1.000      1.0000     0.91107      0.7006     0.80628     0.72551
## Neg Pred Value           1.0000      1.000      1.0000     0.95025      0.9869     0.97676     0.96269
## Prevalence               0.1277      0.118      0.1283     0.12391      0.1272     0.11672     0.13172
## Detection Rate           0.1277      0.118      0.1283     0.07844      0.1163     0.09625     0.09953
## Detection Prevalence     0.1277      0.118      0.1283     0.08609      0.1659     0.11937     0.13719
## Balanced Accuracy        1.0000      1.000      1.0000     0.81215      0.9285     0.89923     0.85613
##                      Class: 80Hz
## Sensitivity              0.60864
## Specificity              0.95367
## Pos Pred Value           0.65559
## Neg Pred Value           0.94387
## Prevalence               0.12656
## Detection Rate           0.07703
## Detection Prevalence     0.11750
## Balanced Accuracy        0.78115

print("Shuffled peaks:")

## [1] "Shuffled peaks:"

# stats
shpk <- caret::confusionMatrix(predictionspk_sh , mu_bs_valid$modulation)
shpk

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction 4Hz 8Hz 16Hz 25Hz 33Hz 40Hz 50Hz 80Hz
##       4Hz  817   0    0    0    0    0    0    0
##       8Hz    0 459   37    0    4    0   88    0
##       16Hz   0  57  412   61  107   59    6   19
##       25Hz   0  25   94  168   60   18   50   40
##       33Hz   0   5   73   59  394   14   29   56
##       40Hz   0  48   66   81   26  357   28   13
##       50Hz   0 134   74  236   33   95  494   57
##       80Hz   0  27   65  188  190  204  148  625
## 
## Overall Statistics
##                                                
##                Accuracy : 0.5822               
##                  95% CI : (0.57, 0.5943)       
##     No Information Rate : 0.1317               
##     P-Value [Acc > NIR] : < 0.00000000000000022
##                                                
##                   Kappa : 0.522                
##                                                
##  Mcnemar's Test P-Value : NA                   
## 
## Statistics by Class:
## 
##                      Class: 4Hz Class: 8Hz Class: 16Hz Class: 25Hz Class: 33Hz Class: 40Hz Class: 50Hz
## Sensitivity              1.0000    0.60795     0.50183     0.21185     0.48403     0.47791     0.58600
## Specificity              1.0000    0.97715     0.94461     0.94881     0.95775     0.95365     0.88681
## Pos Pred Value           1.0000    0.78061     0.57143     0.36923     0.62540     0.57674     0.43989
## Neg Pred Value           1.0000    0.94907     0.92798     0.89487     0.92721     0.93254     0.93386
## Prevalence               0.1277    0.11797     0.12828     0.12391     0.12719     0.11672     0.13172
## Detection Rate           0.1277    0.07172     0.06438     0.02625     0.06156     0.05578     0.07719
## Detection Prevalence     0.1277    0.09187     0.11266     0.07109     0.09844     0.09672     0.17547
## Balanced Accuracy        1.0000    0.79255     0.72322     0.58033     0.72089     0.71578     0.73641
##                      Class: 80Hz
## Sensitivity              0.77160
## Specificity              0.85295
## Pos Pred Value           0.43193
## Neg Pred Value           0.96265
## Prevalence               0.12656
## Detection Rate           0.09766
## Detection Prevalence     0.22609
## Balanced Accuracy        0.81228

validation_sh <- list(shtr, shpk)

How did variable importance metrics compare between the full (unshuffled) version and the two shuffled versions?

importance_sh <- importance(shTrough_rf_classifier)  %>% importance_to_df()
importance_shpk <- importance(shPeak_rf_classifier) %>% importance_to_df()
compare <- rbind(importance_shpk %>% mutate(shuffle = "shuffled peaks"),
                  importance_sh %>% mutate(shuffle = "shuffled troughs"), 
                  importance %>% mutate(shuffle = "original"))

p.imp.sh <- p.importance_shuffled(compare)
p.imp.sh

ggsave(file.path(figures_path,paste0(exp, '_fig3c.png')),p.imp.sh, width = 4, height = 4)

Plot shuffled confusion matrices from training models

Finally, here’s how our shuffled modelled fared at classifying modulation rate classes…

# validation matrices
confmats_sh <- lapply(validation_sh, function(x) cormat_to_df(x$table %>% as.data.frame()) ) %>%
  setNames(c("troughs", "peaks")) %>%
  map_df(., ~as.data.frame(.x), .id="model") %>%
  mutate(model = factor(model, levels = c("troughs", "peaks", "full")))


confmat_sh <- p.mat(confmats_sh, ulim = max(confmats_sh$n)) + 
  facet_grid(cols=vars(model)) 

  print(confmat_sh)

  ggsave(file.path(figures_path,paste0(exp, '_figs6.png')), 
         confmat_sh, width = h, height = w3)

Quantify overlapping calls

overlap_calls <- data[which(data$onset_interval<data$duration),]%>% 
  mutate(overlap = duration-onset_interval)

cat("Of a total of", nrow(data), "calls,", nrow(overlap_calls),
    "(", round(nrow(overlap_calls)/nrow(data)*100,2), "%) overlapped in time.")

## Of a total of 1088994 calls, 31942 ( 2.93 %) overlapped in time.

n_overlap <- overlap_calls %>% group_by(modulation, condition) %>% summarise(n = n())

(n_overlap_conttable <- addmargins(table(overlap_calls$modulation,overlap_calls$condition), c(1,2)))

##       
##        silence steady-state masker random masker   Sum
##   4Hz     3514                 406          1331  5251
##   8Hz     3461                1322           603  5386
##   16Hz    2452                 847           330  3629
##   25Hz    3278                 383           221  3882
##   33Hz    3506                 345           164  4015
##   40Hz    2988                 338           130  3456
##   50Hz    3017                 219            99  3335
##   80Hz    2687                 233            68  2988
##   Sum    24903                4093          2946 31942

Modelling call overlap incidence

Let’s run the overlap data through a negative binomial…

Deviance tables…

(ovlp.anova <- do.call(rbind, lapply(ovlp.models, '[[', 'deviance')) %>% 
   mutate(across(is.numeric, round, 2)))

##            LR Chisq Df Pr(>Chisq)
## condition      5.23  2       0.07
## condition1     7.27  2       0.03
## condition2    12.19  2       0.00
## condition3    15.20  2       0.00
## condition4    25.20  2       0.00
## condition5    31.39  2       0.00
## condition6    34.76  2       0.00
## condition7    35.55  2       0.00

And contrasts…

(contrasts <- do.call(rbind, lapply(ovlp.models, '[[', 'contrasts')) %>%
   mutate(across(is.numeric, round, 2)) %>%
   mutate(modulation = factor(rep(unique(n_overlap$modulation),each=length(unique(n_overlap$condition)))),.before=1) %>%
   remove_rownames() %>%
   dplyr::select(-c("df", "null")) %>%
   dplyr::rename(IRR = ratio))

##    modulation                     masking_condition   IRR    se asymp_lcl asymp_ucl z_ratio p_value
## 1         4Hz       silence / (steady-state masker)  4.58  3.04      0.94     22.45    2.29    0.07
## 2         4Hz               silence / random masker  2.33  1.33      0.60      9.11    1.49    0.41
## 3         4Hz (steady-state masker) / random masker  0.51  0.34      0.10      2.58   -1.00    0.96
## 4         8Hz       silence / (steady-state masker)  1.92  1.21      0.43      8.65    1.04    0.90
## 5         8Hz               silence / random masker  5.36  3.16      1.31     21.96    2.85    0.01
## 6         8Hz (steady-state masker) / random masker  2.79  1.78      0.60     12.88    1.61    0.32
## 7        16Hz       silence / (steady-state masker)  2.21  1.16      0.63      7.77    1.52    0.39
## 8        16Hz               silence / random masker  6.56  3.32      1.95     22.02    3.72    0.00
## 9        16Hz (steady-state masker) / random masker  2.96  1.60      0.81     10.83    2.00    0.14
## 10       25Hz       silence / (steady-state masker)  5.14  3.10      1.21     21.80    2.71    0.02
## 11       25Hz               silence / random masker  9.89  5.80      2.43     40.26    3.91    0.00
## 12       25Hz (steady-state masker) / random masker  1.93  1.27      0.40      9.37    0.99    0.96
## 13       33Hz       silence / (steady-state masker)  5.72  3.07      1.58     20.69    3.24    0.00
## 14       33Hz               silence / random masker 14.70  7.47      4.35     49.66    5.28    0.00
## 15       33Hz (steady-state masker) / random masker  2.57  1.50      0.64     10.41    1.62    0.32
## 16       40Hz       silence / (steady-state masker)  6.48  3.29      1.92     21.87    3.68    0.00
## 17       40Hz               silence / random masker 18.39  9.20      5.55     60.93    5.82    0.00
## 18       40Hz (steady-state masker) / random masker  2.84  1.53      0.78     10.35    1.93    0.16
## 19       50Hz       silence / (steady-state masker) 11.27  6.43      2.88     44.16    4.25    0.00
## 20       50Hz               silence / random masker 27.70 15.51      7.26    105.80    5.93    0.00
## 21       50Hz (steady-state masker) / random masker  2.46  1.45      0.60     10.12    1.52    0.38
## 22       80Hz       silence / (steady-state masker)  7.69  3.88      2.30     25.73    4.04    0.00
## 23       80Hz               silence / random masker 26.34 13.56      7.68     90.35    6.35    0.00
## 24       80Hz (steady-state masker) / random masker  3.43  1.94      0.89     13.24    2.18    0.09

And finally, predict values…

ovlp_preds <- do.call(rbind, lapply(ovlp.models, '[[', 'predictions')) %>%
   mutate(across(is.numeric, round, 2)) %>%
   mutate(modulation = factor(rep(unique(n_overlap$modulation),each=length(unique(n_overlap$condition))))) %>%
   remove_rownames()

ovlp_pred <- p.n_pred(ovlp_preds) %>% recolor()

ovlp_pred

ggsave(file.path(figures_path,paste0(exp, '_figs7b.png') ), ovlp_pred, width = 5, height = 4)

Finally, we can calculate circular stats for overlapping calls, and visualize the angular vectors representing those distributions for each condition:

overlap_circsummary <- overlap_calls %>%
    group_by(modulation, condition) %>%
    summarize(n = n(),
           theta_bar = mean.circular(start_phase_circ, na.rm = T),
           r_bar = rho.circular(start_phase_circ, na.rm = T),
           vm = var.circular(start_phase_circ), 
           rt = rayleigh.test(start_phase_circ, mu=circular(0))$statistic,
           p = rayleigh.test(start_phase_circ, mu=circular(0))$p.value,
           bs_mu_ci_l = mle.vonmises.bootstrap.ci(start_phase_circ, alpha=.05)$mu.ci[1], 
           bs_mu_ci_h = mle.vonmises.bootstrap.ci(start_phase_circ, alpha=.05)$mu.ci[2],) %>%
    mutate(p_adjust = p.adjust(p, method ="bonferroni"), .after = "p") %>%
    mutate(across(where(is.numeric), round, 3)) 

overlap_circsummary %>%
  kbl(escape = FALSE, 
      booktabs = TRUE, 
      format = "html",
      caption = "Circular statistics for overlapping calls") %>%
  kable_styling() %>%
  cat(., file = file.path(tables_path,paste0(exp, "_circ_stats_overlap.html")))

p.ang <- p.angular_vectors(overlap_circsummary %>% confint_wrap())

p.ang

 ggsave(file.path(figures_path,paste0(exp, '_figs7a.png') ), p.ang, width = w, height = 4)

po <- p.circ_density(dat=overlap_calls, circ_dat=overlap_circsummary %>% confint_wrap(), bins = 30) %>% recolor(which = "f")
  print(po)

  ggsave(file.path(figures_path,paste0(exp, '_figs7a-alt2.png') ), po, width = w, height = 5)

Investigate the onset of this adaptation ability

To determine how many cycles of continuous stimulation are needed for the bats’s timing adaptation to stabilize, estimate mean and concentration parameters for successively larger cycle windows:

max_n <- data_init %>% group_by(modulation, condition) %>% summarise(n = max(start_cycle)) %>% 
  ungroup() %>% summarise(min = min(n))

by_cycle_summary <- data.frame()
for (c in seq(0, max_n$min, 10)) {
#for (c in 2:length(seq(0, max_n, 10))) {
  #s <- seq(0, 1000, 10)
  dat_summary <- data_init %>% 
    dplyr::filter(start_cycle <= c) %>%
    #dplyr::filter(start_cycle > s[c-1] & start_cycle <= s[c]) %>% # non-cumulative
    group_by(modulation, condition) %>% # pool over groups
    summarize(n = c,
           theta_bar = mean.circular(start_phase_circ, na.rm = T) %>% as.numeric(),
           r_bar = rho.circular(start_phase_circ, na.rm = T))
  by_cycle_summary <- rbind(by_cycle_summary, dat_summary)
}

References

Some references for circular stats can be found from the NCSS website, as well as the documentation for the Matlab library for circular stats, and articles written on Bayesian modelling of circular data.

There is also the textbook, Statistical analysis of circular data, as well as: - Batschelet, E. (1981). Circular Statistics in Biology. Academic Press: London. - Jammalamadaka, S.R. & Sengupta, A. (2001). Topics in Circular Statistics. World Scientific, River Edge, N.J.